Baixe Matemática Superior para Engenharia - Volume 3 - Ed. 9 -Kreyszig, Erwin e outras Notas de estudo em PDF para Engenharia de Produção, somente na Docsity! Análise Numérica Programa Computacional CAPÍTULO 19 Métodos Numéricos em Geral CAPÍTULO 20 Métodos Numéricos de Álgebra Linear CAPÍTULO 21 Métodos Numéricos para EDOs e EDPs A análise numérica preocupa-se com os métodos numéricos, ou seja, os métodos para resolver problemas em termos de números ou de representações gráficas correspondentes. Ela também inclui a investigação da extensão da aplicabilidade, bem como a precisão e a estabilidade desses métodos. Tarefas usuais em análise numérica são o cálculo de valores de integrais definidas, a obtenção de solução de equações e sistemas lineares, a solução de equações diferenciais ou integrais para as quais inexistem fórmulas de solução e o cálculo de valores de dados experimentais para os quais desejamos obter, por exemplo, um polinômio de aproximação. Os métodos numéricos fornecem, portanto, uma transição de um modelo matemático para um algoritmo, ou seja, uma detalhada receita para se resolver um determinado problema a ser resolvido por programação em seu computador, utilizando um sistema de álgebra computacional (SAC) ou algum outro programa, ou mesmo uma calculadora programável. Neste e nos dois capítulos seguintes, explicamos e ilustramos os métodos numéricos fundamentais mais freqüentemente utilizados em forma algorítmica. O Capítulo 19 trata dos métodos numéricos em geral; o Capítulo 20, da álgebra linear numérica e, em particular, dos métodos empregados na solução de sistemas lineares e de problemas de autovalores de matrizes; e o Capítulo 21 trata dos métodos numéricos aplicados às EDOs e EDPs. Os algoritmos são fornecidos de modo a mostrar da melhor forma possível o funcionamento de cada método. Sugerimos que você também faça uso dos programas comerciais ou de domínio público que estão listados ao final deste texto ou que podem ser obtidos na Internet. Os métodos numéricos vêm cada vez mais ganhando importância em engenharia, mais do que em qualquer outro campo da matemática, em virtude do contínuo desenvolvimento de poderosas técnicas de programa- ção resultantes de uma extensa atividade de pesquisa nessa área, através da invenção de novos métodos, do aperfeiçoamento e adaptação de métodos já existentes e também da redescoberta de outros métodos antigos, que não tinham praticidade na era pré-digital. Um dos principais objetivos dessas atividades é o P A R T E E C A P Í T U L O 1 9 Métodos Numéricos em Geral O primeiro capítulo sobre métodos numéricos começa com uma explicação de alguns conceitos gerais, como ponto de flutuação, erros de arredondamento e erros numéricos em geral e sua propagação. Na Seção 19.2, discutiremos métodos para resolver equações. Os métodos de interpolação, incluindo o dos splines, seguem- se nas Seções 19.3 e 19.4. A última seção (19.5) refere-se à integração e à derivação numéricas. O propósito deste capítulo é duplo. Primeiro, para todas estas tarefas, o aluno deverá familiarizar-se com os métodos numéricos de solução mais básicos (embora não por demais complicados). Eles são indis- pensáveis ao engenheiro porque, para muitos problemas, não existe uma fórmula de solução (imagine, por exemplo, uma integral complicada, ou um polinômio de alto grau ou a interpolação de valores obtidos através de medições). Em outros casos, uma fórmula complicada de solução pode existir, porém não pode ser utilizada na prática. O segundo propósito deste capítulo é fazer com que a aprendizagem dos alunos lhes permita compreender algumas idéias e conceitos básicos que se constituem em importantes métodos completos, como a forma prática de algoritmos, a estimativa de erros e a ordem de convergência. Pré-requisito: Cálculo elementar Referências e Respostas dos Problemas: Parte E do Apêndice 1 e Apêndice 2 19.1 Introdução Os métodos numéricos são utilizados para resolver problemas em computadores e calculadoras através de cálculos numéricos, resultando em uma tabela de números e/ou representações gráficas (figuras). As fases que vão de uma dada situação (em engenharia, economia etc.) até a resposta final são geralmente as seguintes: 1. Modelagem. Estabelecemos um modelo matemático para nosso problema, como uma integral, um sistema de equações, ou uma equação diferencial. 2. Escolha de um método numérico e de parâmetros (p. ex., tamanho do passo de iteração), talvez com uma estimativa preliminar de erro. 3. Programação. Usamos o algoritmo para escrever um programa correspondente em um SAC, como Maple, Mathematica, Matlab ou Mathcad, ou, digamos, em Fortran, C, ou C++, selecionando rotinas compatíveis de um sistema de software conforme o necessário. 4. Realização do cálculo. 5. Interpretação dos resultados em termos físicos ou outros, também decidindo reexecutar o programa caso se necessitem de resultados posteriores. As Etapas 1 e 2 estão relacionadas. Uma leve mudança do modelo pode muitas vezes permitir o uso de um método mais eficiente. Mas, para escolher os métodos, precisamos antes conhecê-los. Os Capítulos 19–21 contêm algoritmos eficientes para as classes mais importantes de problemas que ocorrem com freqüência na prática. Na Etapa 3, o programa consiste em dados fornecidos e em uma seqüência de instruções a serem executadas pelo computador em determinada ordem para produzir a resposta em forma numérica ou gráfica. Para criar uma boa compreensão da natureza do trabalho numérico, continuaremos nesta seção com algumas observações gerais simples. Números em Forma de Ponto Flutuante Sabemos que, na notação decimal, cada número real é representado por uma seqüência finita ou infinita de dígitos decimais. Hoje, a maioria dos computadores tem duas maneiras de representar os números, chamadas de ponto fixo e ponto flutuante. Em um sistema de ponto fixo, todos os números são expressos com um número fixo de Capítulo 19: Métodos Numéricos em Geral 5 decimais à direita da vírgula decimal; por exemplo: números expressos com 3 decimais são 62,358; 0,014; 1,000. Em um texto, escreveríamos, digamos, 3 decimais como 3D. As representações de ponto fixo não são práticas na maioria das computações científicas por causa de seu alcance limitado (explique!) e não nos interessarão. Em um sistema de ponto flutuante, escrevemos, por exemplo: 0,6247 103, 0,1735 1013, 0,2000 101 ou, às vezes, também 6,247 102, 1,735 1014, 2,000 102. Vemos que, neste sistema, o número de dígitos significativos é mantido fixo, ao passo que a vírgula decimal é “flutuante”. Aqui, o dígito significativo de um número c é qualquer dígito dado de c, exceto, possivelmente, para zeros à esquerda do primeiro dígito não-zero; esses zeros servem apenas para fixar a posição da vírgula decimal. (Assim, qualquer outro zero é um dígito significativo de c.) Por exemplo: cada um dos números 1360, 1,360, 0,001360 tem 4 dígitos significativos. Em um texto, indicamos, por exemplo, 4 dígitos significativos por 4S. O uso de expoentes permite-nos representar números muito grandes e números muito pequenos. De fato, teoricamente qualquer número não-zero a pode ser escrito como (1) a m 10n, 0,1 m 1, n inteiro. No computador, m limita-se a k dígitos (p. ex.: k = 8) e n é limitado, dando as representações (apenas para uma quantidade finita de números!) (2) a m 10 n, m 0.d1d2 • • • dk, d1 0. Esses números ā são freqüentemente chamados de números de máquina decimais de k dígitos. Sua parte fra- cionária m (ou m) é chamada de mantissa. Isto nada tem a ver com a “mantissa” que se usa em logaritmos. n é chamado de expoente de ā. Underflow e Overflow. O intervalo de expoentes com que um computador típico pode lidar é muito grande. O ponto flutuante padrão do IEEE (Institute of Electrical and Electronics Engineers) para a precisão simples (o número usual de dígitos nos cálculos) é cerca de –38 < n < 38 (cerca de –125 < n* < 125 para o expoente em representações binárias, ou seja, representações de base 2). [Para a chamada precisão dupla, é cerca de –308 < n < 308 (cerca de –1020 < n* < 1020 para representações binárias)]. Se, em um cálculo ocorrer um número fora desse intervalo faixa, isto é chamado de underflow quando o número é menor e de overflow quando ele é maior. No caso de underflow, o resultado é normalmente arredondado para zero e o cálculo prossegue. O overflow faz o computador parar. Os códigos-padrão (de IMSL, NAG etc.) são escritos de modo a se evitar o overflow. As men- sagens de erro sobre overflow podem então indicar erros de programação (dados incorretos de entrada etc.). Arredondamento Provocam-se erros por causa do corte (= o descarte de todos os números decimais de uma certa casa em diante) ou do arredondamento (aproximação). Esses erros são chamados de erros de arredondamento, independente- mente de terem sido causados pelo corte ou pelo arredondamento propriamente dito. A regra para se arredondar um número para k decimais é a seguinte. (A regra de arredondamento para k dígitos significativos é a mesma, bastando-se substituir a expressão “decimal” por “dígito significativo”.) Regra de Arredondamento. Descarte o (k + 1)-ésimo decimal e todos os decimais subseqüentes. (a) Se a parte assim descartada for inferior à metade da unidade situada na k-ésima casa, deixe o k-ésimo decimal inalterado (“arredondamento para baixo”). (b) Se a parte descartada for maior que a metade da unidade situada na k-ésima casa, aumente de uma unidade a k-ésima casa (“arredondamento para cima”). (c) Se a parte descartada for exatamente igual à meia unidade, arredonde para o decimal par mais próximo. (Exemplo: no arredondamento de 3,45 e 3,55 para uma casa decimal, obtemos 3,4 e 3,6, respectivamente.) O propósito da última parte desta regra é garantir que, nos casos onde a parte descartada vale exatamente a metade de um decimal, os arredondamentos para mais e para menos venham, em média, a ser feitos com apro- ximadamente a mesma freqüência. Se arredondamos 1,2535 para 3, 2, 1 casas decimais, teremos 1,254; 1,25 e 1,3, mas se 1,25 for arredondado para uma casa decimal, sem informações adicionais, teremos 1,2. 6 Parte E • Análise Numérica Simplesmente cortar os números à direita de uma certa casa não é recomendável, pois o erro correspondente pode ser maior do que o ocorrido no arredondamento, além de ser sistemático. (Não obstante, alguns computadores utilizam tal procedimento, pelo fato de ele ser mais simples e mais rápido. Por outro lado, alguns computadores e calculadoras aumentam a precisão de seus resultados fazendo cálculos intermediários que utilizam um ou mais dígitos extras, chamados de dígitos de guarda.) Erro de Arredondamento. Em (2), chamemos de a = fl(a) a aproximação de ponto flutuante que o computador faz de a em (1) utilizando o arredondamento, onde fl sugere flutuante. Então a aplicação da regra do arredondamento fornece (eliminando-se os expoentes) m m 12 10 k. Como m 0,1, isto implica que (sendo a 0) (3) a a a m m m 1 2 101k. O lado direito u 1 2 101k desta equação é chamado de unidade de arredondamento. Se escrevermos a = a(1 + d), teremos, usando álgebra, que (a a)/a d, logo, d u por (3). Isto mostra que a unidade de arre- dondamento u é um limite de erro do arredondamento. Os erros de arredondamento podem arruinar completamente um cálculo, mesmo quando ele for pequeno. Em geral, esses erros vão se tornando cada vez mais perigosos à medida que se aumenta o número de operações arit- méticas a serem executadas (com esse número chegando às vezes a vários milhões!). Portanto, é importante que analisemos os programas computacionais com o propósito de descobrir erros esperados de arredondamento e de encontrar para os cálculos um arranjo que permita que os erros de arredondamento sejam os menores possíveis. Tampouco a aritmética dos computadores é exata e ela provoca outros erros; contudo, esses erros não serão relevantes à nossa discussão. Precisão nas Tabelas. Embora os programas disponíveis hoje em dia tenham tornado supérfluas diversas tabelas de valores de funções, algumas tabelas (de funções superiores, de coeficientes de fórmulas de integração etc.) ainda continuarão sendo de uso ocasional. Se uma tabela apresenta k dígitos significativos, presume-se, por convenção, que qualquer valor a na tabela se desvia do valor exato a por, no máximo, 1 2 unidade do k-ésimo dígito. Algoritmo. Estabilidade Os métodos numéricos podem ser formulados como algoritmos. Um algoritmo é um processamento passo a passo que apresenta um método numérico sob uma forma (um “pseudocódigo”) compreensível aos seres humanos. (Vire algumas páginas para ver a aparência que tem um algoritmo.) Os algoritmos são então utilizados para se escreverem programas em uma linguagem de programação que seja compreendida pelo computador, de modo que este possa executar o método numérico. As seções a seguir apresentam alguns algoritmos importantes. Para as tarefas de rotina, pode ocorrer que seu sistema de álgebra computacional ou algum outro aplicativo contenha programas que você possa utilizar ou incluir como partes de programas que você mesmo criar. Estabilidade. Para ser útil, um algoritmo deve ser estável; ou seja, pequenas alterações nos dados iniciais devem causar apenas pequenas alterações nos resultados finais. Entretanto, se pequenas alterações nos dados iniciais provocarem grandes alterações nos resultados finais, chamaremos o algoritmo de instável. É preciso distinguir essa “instabilidade numérica”, que na maioria dos casos pode ser evitada com a escolha de um algoritmo melhor, da “instabilidade matemática” de um problema, a qual é também chamada de “mau condicionamento”, um conceito que discutiremos na próxima seção. Alguns algoritmos são estáveis somente para certos dados iniciais, de modo que é preciso que se tome cuidado em casos assim. Erros de Resultados Numéricos Os resultados finais de cálculos de quantidades desconhecidas geralmente são aproximações; ou seja, não são exatos, mas envolvem erros, os quais podem resultar de uma combinação dos seguintes efeitos: os erros de arre- dondamento, resultantes dos procedimentos de arredondamento, conforme já discutimos; os erros experimentais, encontrados nos dados fornecidos (provavelmente surgindo das medições); os erros de truncamento, resultantes dos procedimentos de truncamento (interrupção prematura), por exemplo, quando substituímos uma série de Taylor pela soma de seus primeiros termos. Esses erros dependem do método computacional empregado e devem ser tratados individualmente, dependendo de cada método. [Às vezes, o termo “truncar” é usado como um sinônimo de “cortar” (ver o que já dissemos deste último), uma terminologia que não recomendamos.] Capítulo 19: Métodos Numéricos em Geral 9 dondar os números arredondados e não os números original- mente dados. 4. Repita o Problema 3 com números de sua própria escolha e que forneçam resultados com diferenças ainda mais drásticas. Como pode você evitar essas dificuldades? 5. O quociente do Problema 3 está na forma a/(b – c). Escreva-o como a(b + c)/(b2 – c2). Calcule-o primeiro com 5S, depois arredondando o numerador 12,996 e o denominador 2,28 paulatinamente como no Problema 3. Compare e comente. 6. (Equação quadrática) Resolva x2 – 20x + 1 = 0 por (6) e por (7), usando 6S no cálculo. Compare e comente. 7. Faça os cálculos do Problema 6 com 4S e 2S. 8. Resolva: x2 – 100x + 2 = 0 por (6) e (7) com 5S e compare. 9. Calcule: 1/e = 0,367879 (6S) usando as somas parciais de 5 a 10 termos da série de Maclaurin: (a) de e–x com x = 1, (b) de ex com x = 1 e calculando depois o recíproco. Qual é o resultado mais preciso? 10. A adição com um número fixo de dígitos significativos depende da ordem em que os números são somados. Ilustre esse fato com um exemplo. Descubra uma regra empírica para a melhor ordem. 11. As aproximações de = 3,141 592 653 589 79 • • • são 22/7 e 355/113. Determine com 3 dígitos significativos os erros corres- pondentes e os erros relativos. 12. Compute p usando a aproximação de Machin 16 arctan (1/5) – 4 arctan (1/239) para 10S (que é um resultado correto). (Em 1986, D. H. Barley calculou quase 30 milhões de decimais de p em um CRAY-2 em menos de 30 horas. E a corrida para encontrar cada vez mais decimais ainda continua.) 13. (Arredondamento e soma). Consideremos que a1, • • •, an sejam números com os aj arredondados corretamente até Dj decimais. Ao calcularmos a soma a1 + • • • + an, retendo D = min Dj deci- mais, qual opção é essencial: primeiro somar e depois arredondar o resultado, ou primeiro arredondar cada número para D decimais e depois somar? 14. (Teoremas sobre erros) Prove o Teorema 1(a) para a adição. 15. Prove o Teorema 1(b) para a divisão. 16. No Exemplo 1, mostre que o valor absoluto do erro de x2 = 2,000/39,95 = 0,05006 é inferior a 0,00001. 17. Overflow e underflow podem às vezes ser evitados por meio de simples mudanças em uma fórmula. Explique isso em termos de x2 y2 x1 (y/x)2 com x2 y2 e um x tão grande que x2 causaria overflow. Crie seus próprios exemplos. 18. (Forma aninhada) Calcule o valor de: ƒ(x) x3 7,5x2 11,2x 2,8 ((x 7,5)x 11,2)x 2,8 em x = 3,94 usando uma aritmética de 3S e arredondando, nas duas formas dadas. A última forma, chamada de aninhada, é geralmente preferível, visto que ela minimiza o número de operações e, con- seqüentemente, o efeito do arredondamento. 19. EXPERIMENTO DE ÁLGEBRA COMPUTACIONAL. Corte e Arredondamento. (a) Considere x = 4/7 e y = 1/3. Encontre os erros ecorte, earred e os erros relativos er,cor, er,arred ocorridos devido aos cortes e aos arredondamentos para 5S em x y, x y, xy, x /y. Experimente com outras frações de sua escolha. (b) Represente em um mesmo gráfico ecorte e earred (até 5S) de k/21 como funções de k = 1, 2, • • • 21. A partir do gráfico, que valor médio se pode perceber para ecorte? E para earred? Experimente com outros números inteiros que forneçam gráficos similares. Idem para os que forneçam diferentes tipos de gráficos. Você poderia caracterizar os diferentes tipos em termos de fatores primos? (c) Como muda a situação em (b) se você escolher 4S em vez de 5S? (d) Escreva programas para os trabalhos em (a)–(c). 20. PROJETO ESCRITO. Métodos numéricos. Usando suas pró- prias palavras, escreva, sobre o papel geral dos métodos numéricos na matemática aplicada, o motivo pelo qual eles são importantes, onde e quando precisam ou podem ser usados, e como o uso do computador os influencia na engenharia e em outras áreas. 19.2 Solução de Equações por Iteração Daqui em diante, cada seção será dedicada a algum tipo básico de problema e aos métodos correspondentes de solução. Começaremos com os métodos que permitem encontrar soluções de uma equação única (1) f(x) = 0 onde f é uma função dada. Para esta tarefa, praticamente não existem fórmulas (exceto em alguns casos simples), de modo que dependemos quase inteiramente dos algoritmos numéricos. Uma solução de (1) é um número x = s, tal que f(s) = 0. Aqui, s sugere “solução”, embora também utilizaremos outras letras. Exemplos de equações assim são x3 + x = 1, sen x = 0,5x, tan x = x, cosh x = sec x, cosh x cos x = –1, que podem ser todas escritas na forma (1). O primeiro desses exemplos refere-se a uma equação algébrica, visto que a função f correspondente é um polinômio e, neste caso, as soluções são também chamadas de raízes da equação. Os demais exemplos são equações transcendentais, pois envolvem funções transcendentais. Resolver as equa- ções (1) é uma tarefa de fundamental importância, visto que elas são abundantes nas aplicações em engenharia: algumas ocorrem nos Capítulos 2, 4, 8 (equações características), 6 (frações parciais), 12 (autovalores, zeros das funções de Bessel) e 16 (integração), embora haja também muitas e muitas mais. Para resolvermos (1) quando não há nenhuma fórmula capaz de fornecer a solução exata, podemos usar um método de aproximação, em particular, um método de iteração, ou seja, um método em que partimos de um palpite inicial x0 (que pode ser ruim) e passo a passo calculamos aproximações x1, x2, • • • (que, de modo geral, 10 Parte E • Análise Numérica são progressivamente melhores) de uma solução desconhecida de (1). Discutiremos aqui três desses métodos que são particularmente importantes do ponto de vista prático e mencionaremos dois outros no conjunto de problemas propostos. Esses métodos e os princípios a eles subjacentes são fundamentais para a compreensão dos diversos métodos empregados nos aplicativos computacionais. De modo geral, os métodos iterativos são fáceis de programar porque as operações computacionais são as mesmas em cada etapa — apenas os dados mudam de uma etapa para outra — e, o que é mais importante, se em um caso concreto um método converge, ele geralmente é estável (veja a Seção 19.1). Iteração de Ponto Fixo para a Resolução de Equações f(x) = 0 O uso que agora estamos fazendo da expressão “ponto fixo” não tem absolutamente nada a ver com o sentido adotado na última seção. De um modo qualquer, transformamos (1) algebricamente na forma (2) x = g(x). Então, escolhemos um x0 e calculamos x1 = g(x0), x2 = g(x1) e, de maneira geral, (3) xn 1 g(xn) (n 0, 1, • • •), Chamamos uma solução de (2) de um ponto fixo de g, o que justifica o nome do método. Esta é também uma solução de (1), visto que, de x = g(x), podemos retornar à forma original f(x) = 0. A partir de (1), podemos obter diferentes formas de (2). O comportamento das seqüências iterativas correspondentes x0, x1,• • • pode diferir, em particular, com respeito à sua velocidade de convergência. De fato, algumas delas podem simplesmente não convergir. Ilustremos esses fatos com um exemplo simples. EXEMPLO 1 Um Processo Iterativo (Iteração de Ponto Fixo) Monte um processo de iteração para a equação f(x) = x2 – 3x + 1 = 0. Como conhecemos as soluções x 1,5 1,25, portanto, 2,618 034 e 0,381 966, podemos observar o comportamento do erro à medida que a iteração prossegue. Solução. A equação pode ser escrita como (4a) x g1(x) 1 3(x 2 1), portanto, xn 1 1 3(xn 2 1), Se escolhermos x0 = 1, obteremos a seqüência (Fig. 423a; calculada com 6S e depois arredondada) x0 1,000, x1 0,667, x2 0,481, x3 0,411, x4 0,390, • • • que parece se aproximar da solução menor. Se escolhermos x0 = 2, a situação será semelhante. Se escolhermos x0 = 3, obteremos a seqüência (Fig. 423a; parte superior) x0 3,000, x1 3,333, x2 4,037, x3 5,766, x4 11,415, • • • que diverge. 0 0 5 5x 0 0 5 5x g 1 (x) g 2 (x) (a) (b) Fig. 423. Exemplo 1, iterações (4a) e (4b) Capítulo 19: Métodos Numéricos em Geral 11 Nossa equação também pode ser escrita (dividindo-se por x) como (4b) x g2(x) 3 1 x , portanto, xn 1 3 1 xn , e, se escolhermos x0 = 1, obteremos a seqüência (Fig. 423b) x0 1,000, x1 2,000, x2 2,500, x3 2,600, x4 2,615, • • • que parece se aproximar da solução maior. De modo semelhante, se escolhermos x0 = 3, obteremos a seqüência (Fig. 423b) x0 3,000, x1 2,667, x2 2,625, x3 2,619, x4 2,618, • • • . Nossos gráficos mostram o seguinte: na parte inferior da Fig. 423a, a inclinação de g1(x) é menor que a inclinação de y = x, a qual vale 1; portanto, g 1(x) 1 e constatamos que há convergência. Na parte superior, g1(x) é mais inclinada (g 1(x) 1) e ocorre divergência. Na Fig. 423b, a inclinação de g2(x) é menor próxima ao ponto de intersecção (x = 2,618, ponto fixo de g2, solução de f(x) = 0) e ambas as seqüências parecem convergir. De tudo isso, concluímos que a convergência parece depender do fato de a curva de g(x) ser menos inclinada que a reta y = x na vizinhança de uma solução, e veremos agora que esta condição g (x) < 1 (= inclinação de y = x) é suficiente para a convergência. Dizemos que um processo de iteração definido por (3) é convergente para um x0 se a seqüência correspondente x0, x1, • • • for convergente. Uma condição suficiente para a convergência é dada no seguinte teorema, que possui diversas aplicações práticas. TEOREMA 1 Convergência da Iteração de Ponto Fixo Consideremos que x = s seja uma solução de x = g(x) e suponhamos que g tenha uma derivada contínua em algum intervalo J contendo s. Então, se g (x) K 1 em J, o processo de iteração definido por (3) converge para qualquer x0 em J e o limite da seqüência {xn} é s. PROVA Pelo teorema do valor médio no cálculo diferencial, existe um t entre x e s tal que g(x) g(s) g (t) (x s) (x em J). Como g(s) = s e x1 = g(x0), x2 = g(x1), • • •, obtemos disso e da condição para g (x) no teorema: xn s g(xn1) g(s) g (t)xn1 s K xn1 s. Aplicando esta desigualdade n vezes, para n, n – 1, • • •, 1 obtemos xn s K xn1 s K2xn2 s • • • Knx0 s. Como K 1, temos Kn → 0; logo, xn s → 0 quando n → . Mencionemos que uma função g que satisfaz à condição no Teorema 1 é chamada de contração, visto que g(x) g(v) Kx v, onde K 1. Além disso, K fornece informações sobre a velocidade da convergência. Por exem- plo, se K = 0,5, então em apenas 7 iterações a precisão aumenta pelo menos em 2 dígitos, pois 0,57 < 0,01. EXEMPLO 2 Um Processo Iterativo. Ilustração do Teorema 1 Encontre uma solução de f(x) = x3 + x = 0 por iteração. Solução. Um esboço mostra que há uma solução próxima de x = 1. Podemos escrever a equação como (x2 + 1)x = 1 ou x g1(x) 1 1 x2 , de modo que xn 1 1 1 xn2 . Também g 1(x) 2x (1 x2)2 1 para qualquer x, pois 4x2/(1 x2)4 4x2/(1 4x2 • • •) 1, de modo que, pelo Teorema 1, temos a convergência para qualquer x0. Escolhendo x0 = 1, obtemos (Fig. 424) x1 0,500, x2 0,800, x3 0,610, x4 0,729, x5 0,653, x6 0,701, • • • . A solução exata para 6S é s = 0,682 328. A equação dada também pode ser escrita como x g2(x) 1 x 3. Então, g 2(x) 3x2 e isso é maior do que 1 próximo da solução, de modo que não podemos aplicar o Teorema 1 e assegurar a convergência. Experimente com x0 = 1; x0 = 0,5; x0 = 2 e veja o que acontece. Este exemplo mostra que a transformação de uma dada f(x) = 0 para a forma x = g(x) com g satisfazendo a g (x) K 1 pode precisar de alguma experimentação. 14 Parte E • Análise Numérica EXEMPLO 5 Método de Newton Aplicado a uma Equação Algébrica Aplique o método de Newton à equação f(x) = x3 + x – 1 = 0. Solução. De (5), temos xn 1 xn xn 3 xn 1 3xn 2 1 2xn 3 1 3xn 2 1 , Começando de x0 = 1, obtemos x1 0,750 000, x2 0,686 047, x3 0,682 340, x4 0,682 328, • • • onde x4 tem o erro –1 10–6. Uma comparação com o Exemplo 2 mostra que a convergência que temos agora é muito mais rápida. Isso pode motivar o conceito de ordem de um processo iterativo, que discutiremos a seguir. Ordem de um Método Iterativo. Rapidez da Convergência A qualidade de um método iterativo pode ser caracterizada pela rapidez da convergência, como se segue. Consideremos que xn+1 = g(xn) defina um método de iteração e que xn seja uma aproximação da solução s de x = g(x). Então, xn = s – en, onde en é o erro de xn. Suponha que g seja derivável várias vezes, de modo que a fórmula de Taylor forneça (6) xn 1 g(xn) g(s) g (s)(xn s) 12g (s)(xn s) 2 • • • g(s) g (s)en 12g (s)en 2 • • • . O expoente de en no primeiro termo após g(s) que não desaparece é chamado de ordem do processo iterativo definido por g. A ordem mede a velocidade da convergência. Para ver isso, subtraia g(s) = s de ambos os lados de (6). Então, no lado direito, faça xn+1 – s = –en+1, onde en+1 é o erro de xn+1, e, no lado direito, a expressão restante é aproximadamente igual ao seu primeiro termo não-nulo, visto que en é pequeno no caso de haver convergência. Portanto, (7) (a) en 1 g (s)en (b) en 1 12g (s)en 2 no caso de primeira ordem, no caso de segunda ordem, etc. Assim, se en = 10–k em alguma iteração, então para a segunda ordem, en 1 c (10k)2 c 102k, de modo que o número de dígitos significativos é aproximadamente duplicado em cada iteração. Convergência do Método de Newton No método de Newton, g(x) x ƒ(x) /ƒ (x). Por derivação, (8) g (x) 1 ƒ (x)2 ƒ(x)ƒ (x) ƒ (x)2 ƒ(x)ƒ (x) ƒ (x)2 . Como f(s) = 0, isto mostra que também g (s) = 0. Logo, o método de Newton é pelo menos de segunda ordem. Se derivarmos novamente e fizermos x = s, descobrimos que (8*) g (s) ƒ (s) ƒ (s) que não será zero em geral. Isto prova o TEOREMA 2 Convergência de Segunda Ordem do Método de Newton Se f(x) é três vezes derivável e f e f não são nulas numa solução s de f(x) = 0, então, para x0 suficiente- mente próximo de s, o método de Newton é de segunda ordem. Comentários. Para o método de Newton, (7b) torna-se, por (8*), Capítulo 19: Métodos Numéricos em Geral 15 (9) en 1 ƒ (s) 2ƒ (s) en2. Para que o método indicado no Teorema 2 apresente uma rápida convergência, é importante que s seja um zero simples de f(x0) (assim, f (s) 0) e que x0 esteja próximo de s, visto que, na fórmula de Taylor, tomamos apenas o termo linear [veja (5*)], supondo que o termo quadrático seja desprezivelmente pequeno. (Com um valor ruim para x0, este método pode até mesmo divergir!) EXEMPLO 6 Estimativa a Priori do Erro para o Número de Etapas da Iteração de Newton Use x0 = 2 e x1 = 1,901 no Exemplo 4 para estimar quantas etapas de iteração são necessárias para se produzir uma solução com uma pre- cisão de 5D. Esta é uma estimativa a priori ou uma estimativa prévia, porque só podemos calculá-la depois de uma iteração e antes das demais iterações. Solução. Temos f(x) = x – 2 sen x = 0. Derivando, obtemos ƒ (s) 2ƒ (s) ƒ (x1) 2ƒ (x1) 2 sen x1 2(1 2 cos x1) 0,57. Logo, (9) fornece en 1 0,57en2 0,57(0,57e2n1)2 0,573e4n1 • • • 0,57Me0M 1 5 106 onde M 2n 2n1 • • • 2 1 2n 1 1. Mostramos a seguir que e0 ≈ –0,11. Conseqüentemente, nossa condição passa a ser 0,57M0,11M 1 5 106. Logo, segundo essa estimativa grosseira, n = 2 é o menor n possível, o que está em boa concordância com o Exemplo 4. O valor e0 ≈ –0,11 é obtido de e1 e0 (e1 s) (e0 s) x1 x0 0,10, donde e1 e0 0,10 0,57e02 ou 0,57e02 e0 0,10 0, que dá e0 0,11. Dificuldades no Método de Newton. As dificuldades podem surgir se f (x) for muito pequeno próximo a uma solução s de f(x) = 0, por exemplo, se s for um zero de f(x) de segunda ordem (ou de ordem mais elevada) (de modo que o método de Newton converge apenas linearmente, como se mostra aplicando-se a regra de l’Hôpital (8)). Geometricamente, um pequeno valor de f (x) significa que a tangente a f(x) próximo a s quase coincide com o eixo x (de modo que se possa precisar de uma precisão dupla para se obterem f(x) e f (x) com precisão suficiente). Então, para valores x s afastados de s, ainda podemos ter pequenos valores de função R(s) ƒ(s). Neste caso, dizemos que a equação f(x) = 0 é mal-condicionada. R(s) é chamado de resíduo de f(x) = 0 em s. Assim, um resíduo pequeno garante um erro pequeno de s somente se a equação não for mal-condicionada. EXEMPLO 7 Uma Equação Mal-condicionada ƒ(x) x5 104x 0 é mal-condicionada. x = 0 é uma solução. f (0) = 10–4 é pequeno. Em s = 0,1, o resíduo f(0,1) = 2 10–5 é pequeno, porém o erro –0,1 é maior por um fator de 5000 em módulo; conceba você mesmo um exemplo ainda mais drástico. Método da Secante para a Resolução de f(x) = 0 O método de Newton é muito poderoso, mas tem a desvantagem de que a derivada f pode às vezes ser uma expressão muito mais difícil do que a própria função f, tornando o cálculo do seu valor computacionalmente exigente. Essa situação sugere a idéia de se substituir a derivada pelo quociente de diferença ƒ (xn) ƒ(xn) ƒ(xn1) xn xn1 . Então, ao invés de (5), temos a fórmula do popular método da secante y x P n –1 x n –1 y = f(x) Secante P n x n s x n +1 Fig. 426. Método da secante 16 Parte E • Análise Numérica (10) xn 1 xn ƒ(xn) xn xn1 ƒ(xn) ƒ(xn1) . Geometricamente, na Fig. 426, obtemos uma intersecção entre o eixo x em xn+1 e a secante a f(x) passando pelos pontos Pn–1 e Pn. Precisamos de valores iniciais, x0 e x1. Dessa forma, evitamos o cálculo de derivadas. Pode-se mostrar que a convergência é superlinear (ou seja, mais rápida do que a linear, en 1 const en1,62; veja [E5] no Apêndice 1]), sendo quase quadrática, como o método de Newton. O algoritmo é semelhante ao do método de Newton, como o aluno pode mostrar. CUIDADO! Não é bom escrever (10) como xn 1 xn1ƒ(xn) xnƒ(xn1) ƒ(xn) ƒ(xn1) , pois isso pode causar uma perda de dígitos significativa se xn e xn-1 forem aproximadamente iguais entre si. (Você consegue ver isso a partir da fórmula?) EXEMPLO 8 Método da Secante Encontre a solução positiva de f(x) = x – 2 sen x = 0 pelo método da secante, começando de x0 = 2, x1 = 1,9. Solução. Aqui, (10) é xn 1 xn (xn 2 sen xn)(xn xn1) xn xn1 2(sen xn1 sen xn) xn Nn Dn . Os valores numéricos são: n xn1 xn Nn Dn xn 1 xn 1 2,000 000 1,900 000 0,000 740 0,174 005 0,004 253 2 1,900 000 1,895 747 0,000 002 0,006 986 0,000 252 3 1,895 747 1,895 494 0 0 x3 = 1,895 494 é exato em 6D. Veja o Exemplo 4. Resumo dos Métodos. Os métodos para calcular soluções s de f(x) = 0 com funções f(x) contínuas (ou deriváveis) dadas começam com uma aproximação inicial x0 de s e geram uma seqüência x1, x2, • • • por iteração. Os métodos de ponto fixo resolvem f(x) = 0 escrita como x = g(x), de modo que s é um ponto fixo de g, ou seja, s = g(s). Para g(x) = x – f(x)/f (x), este é o método de Newton, o qual, para um bom valor de x0 e zeros simples da função, converge quadraticamente (e converge linearmente quando há múltiplos zeros). O método da secante decorre do de Newton substituindo-se f (x) por um quociente de diferenças. O método da bissecção e o método de falsa posição, apresentados em Problemas Propostos 19.2, sempre convergem, embora freqüentemente devagar. P R O B L E M A S P R O P O S T O S 1 9 . 2 1–7 ITERAÇÃO DE PONTO FIXO Aplique a iteração de ponto fixo e responda às perguntas relacionadas quando indicado. Mostre os detalhes do que fizer. 1. x = 1,4 sen x, x0 = 1,4 2. Faça as iterações indicadas no final do Exemplo 2. Esboce uma figura semelhante à Fig. 424. 3. Por que obtemos uma seqüência monótona no Exemplo 1, mas não no Exemplo 2? 4. f = x4 – x + 0,2 = 0, raiz próxima de 1, x0 = 1 5. f como no Problema 4, raiz próxima de 0, x0 = 0 6. Encontre a menor solução positiva de sen x = e–0,5x, x0 = 1. 7. (Funções de Bessel, batida do tambor) Uma soma parcial da série de Maclaurin de J0(x) (Seção 5.5) é ƒ(x) 1 1 4x 2 1 64x 4 12304x 6. Conclua, a partir de um esboço, que, próximo de x = 2, f(x) = 0. Escreva f(x) = 0 como x = g(x) (dividindo f(x) por 14x e passando o termo em x resultante para o outro lado). Encontre o zero. (Veja na Seção 12.9 a importância desses zeros.) 8. PROJETO DE ÁLGEBRA COMPUTACIONAL. Iteração de Ponto Fixo. (a) Existência. Prove que se g é contínua em um intervalo fechado I e seu intervalo de definição se situa em I, então a equação x = g(x) tem pelo menos uma solução em I. Ilustre o fato de que ela pode ter mais de uma solução em I. b) Convergência. Considere que ƒ(x) x3 2x2 3x 4 0. Escreva isto como x = g(x), escolhendo g como (1) (x3 ƒ)1/3, (2) (x2 12ƒ) 1/2, (3) x 13ƒ, (4) x(1 1 4ƒ), (5) (x 3 ƒ)/x2, (6) (2x2 ƒ)/2x, (7) x ƒ/ƒ e, em cada caso, usando x0 = 1,5. Investigue a convergência e a divergência, bem como o número de iterações necessárias para se alcançarem valores exatos em 6S para raízes. Capítulo 19: Métodos Numéricos em Geral 19 L0(x) x 9,5 0,5 2,0(x 9,5), L0(9,2) 2,0(0,3) 0,6 L1(x) x 9,0 0,5 2,0(x 9,0), L1(9,2) 2 0,2 0,4 (veja a Fig. 429) e obtemos a resposta ln 9,2 p1(9,2) L0(9,2)ƒ0 L1(9,2)ƒ1 0,6 2,1972 0,4 2,2513 2,2188. O erro é e a a 2,2192 2,2188 0,0004. Logo, a interpolação linear não é aqui suficiente para obtermos uma precisão de 4D; ela teria sido suficiente para uma precisão de 3D. 9 9,59,2 10 11 x 0 1 y L0 L1 Fig. 429. L0 e L1 do Exemplo 1 A interpolação quadrática é a calculada para os pontos dados: (x0, f0), (x1, f1), (x2, f2), usando-se um polinômio de segundo grau p2(x), o qual, pela idéia de Lagrange, corresponde a (3a) p2(x) L0(x)ƒ0 L1(x)ƒ1 L2(x)ƒ2 com L0(x0) 1, L1(x1) 1, L2(x2) 1, e L0(x1) L0(x2) 0 etc. Estabelecemos que L0(x) l0(x) l0(x0) (x x1)(x x2) (x0 x1)(x0 x2) (3b) L1(x) l1(x) l1(x1) (x x0)(x x2) (x1 x0)(x1 x2) L2(x) l2(x) l2(x2) (x x0)(x x1) (x2 x0)(x2 x1) . Como obtivemos este resultado? Ora, o numerador faz Lk(xj) = 0 se j k. E o denominador faz Lk(xk) = 1, porque é igual ao numerador em x = xk. EXEMPLO 2 Interpolação Quadrática de Lagrange Calcule ln 9,2 por (3), usando os dados fornecidos no Exemplo 1 e um terceiro valor adicional ln 11,0 = 2,3979. Solução. Em (3), L0(x) (x 9,5)(x 11,0) (9,0 9,5)(9,0 11,0) x2 20,5x 104,5, L0(9,2) 0,5400, L1(x) (x 9,0)(x 11,0) (9,5 9,0)(9,5 11,0) 1 0,75 (x2 20x 99), L1(9,2) 0,4800, L2(x) (x 9,0)(x 9,5) (11,0 9,0)(11,0 9,5) 1 3 (x2 18,5x 85,5), L2(9,2) 0,0200, (veja a Fig. 430), de modo que (3a) fornece, com um valor exato em 4D, ln 9,2 p2(9,2) 0,5400 2,1972 0,4800 2,2513 0,0200 2,3979 2,2192. 0 1 9 9,5 10 11 x y L0 L2 L1 Fig. 430. L0, L1, L2 do Exemplo 2 20 Parte E • Análise Numérica Polinômio Geral de Interpolação de Lagrange. Para n geral, obtemos (4a) ƒ(x) pn(x) n k0 Lk(x)ƒk n k0 lk(x) lk(xk) ƒk onde Lk(xk) = 1, Lk é 0 nos outros nós e os Lk são independentes da função f a ser interpolada. Obtemos (4a) se fizermos l0(x) (x x1)(x x2) • • • (x xn), (4b) lk(x) (x x0) • • • (x xk1)(x xk 1) • • • (x xn), 0 k n, ln(x) (x x0)(x x1) • • • (x xn1). Podemos ver facilmente que pn(xk) = fk. De fato, inspecionando (4b), vemos que lk(xj) = 0 se j k, de modo que, para x = xk, a soma em (4a) se reduz ao único termo (lk(xk)/lk(xk))fk = fk. Estimativa de Erro. Se a própria função f for um polinômio de grau n (ou menor), ela deve coincidir com pn, porque os n + 1 dados (x0, f0),• • •, (xn, fn) determinam unicamente um polinômio, de modo que o erro é zero. Ora, a função especial f tem sua (n + 1)-ésima derivada identicamente nula. Isto torna plausível o fato de que, para uma função geral f, sua (n + 1)-ésima derivada f(n+1) deve medir o erro en(x) ƒ(x) pn(x). Pode-se mostrar que isso é verdadeiro se f(n+1) existir e for contínua. Então, com um valor adequado para t situado entre x0 e xn (ou entre x0, xn e x se extrapolarmos), (5) en(x) ƒ(x) pn(x) (x x0)(x x1) • • • (x xn) ƒ(n 1)(t) (n 1)! . Portanto, e(x) é 0 nos nós e pequeno próximo a eles, por causa da continuidade. O produto (x – x0)• • •(x – xn) é grande para valores de x distantes dos nós. Isto torna a extrapolação arriscada. E a interpolação em um valor x será melhor se escolhermos os nós em ambos os lados desse x. Além disso, obtemos os limites de erro tomando o menor e o maior valores de f(n+1)(t) em (5) no intervalo x0 t xn (ou no intervalo que também contém x se extrapolarmos). O ponto mais importante é que, como pn é único, conforme mostramos, temos o TEOREMA 1 Erro de Interpolação A fórmula (5) fornece o erro para qualquer método de interpolação polinomial se f(x) possui uma (n + 1)-ésima derivada contínua. Estimativa prática de erro. Se a derivada em (5) for difícil ou impossível de ser obtida, aplique o Princípio do Erro (Seção 19.1), ou seja, tome outro nó e o polinômio de Lagrange pn+1(x), e considere pn+1(x) – pn(x) como uma estimativa de erro (grosseira) de pn(x). EXEMPLO 3 Estimativa de Erro (5) da Interpolação Linear. Prejuízo por Arredondamento. Princípio do Erro Estime o erro no Exemplo 1, primeiro diretamente por (5) e depois pelo Princípio do Erro (Seção 19.1). Solução. (A) Estimativa por (5). Temos n 1, ƒ(t) ln t, ƒ (t) 1/t, ƒ (t) 1/t2. Logo, e1(x) (x 9,0)(x 9,5) (1) 2t2 , portanto, e1(9,2) 0,03 t2 . t = 9,0 fornece o valor máximo 0,03/92 = 0,00037 e t = 9,5 fornece o valor mínimo 0,03/9,52 = 0,00033, de modo que obtemos 0,00033 e1(9,2) 0,00037, ou melhor, 0,00038, visto que 0,3/81 = 0,003 703 • • • . Entretanto, isto está em discordância com o Exemplo 1, para o qual o erro vale 0,0004, de modo que podemos aprender alguma coisa! Para este caso, a repetição do cálculo com 5D em vez de 4D fornece ln 9,2 p1(9,2) 0,6 2,19722 0,4 2,25129 2,21885 com um erro real e 2,21920 2,21885 0,00035 situado bem próximo do meio entre nossos dois limites de erro. Isto mostra que a discrepância (0,0004 versus 0,00035) foi causada pelo arredondamento, algo que não foi levado em conta em (5). Capítulo 19: Métodos Numéricos em Geral 21 (B) Estimativa pelo Princípio do Erro. Calculamos p1(9,2) = 2,21885 como antes e depois p2(9,2) como no Exemplo 2, porém com 5D, obtendo p2(9,2) 0,54 2,19722 0,48 2,25129 0,02 2,39790 2,21916. A diferença p2(9,2) – p1(9,2) = 0,00031 é o erro aproximado de p1(9,2) que queríamos obter; esta é uma aproximação do verdadeiro erro 0,00035 dado anteriormente. Interpolação das Diferenças Divididas de Newton Para os dados fornecidos (x0, f0),• • •, (xn, fn), o polinômio de interpolação pn(x) que satisfaz a (1) é único, con- forme já demonstramos. Porém, para diferentes propósitos, podemos usar pn(x) expresso em formas diferentes. A forma de Lagrange, recém-discutida, é útil para obter fórmulas de integração (Seção 19.5) e de derivação numéricas (fórmulas de aproximação para derivadas). De maior importância prática são as formas de Newton para pn(x), que também usaremos para resolver EDOs (na Seção 21.2). Elas envolvem menos operações aritméticas que a forma de Lagrange. Além disso, freqüentemente ocorre que temos de aumentar o grau n para alcançarmos uma precisão exigida. Nesses casos, com as formas de Newton, podemos usar todo o trabalho que já desenvolvemos e apenas acrescentar um outro termo, uma possibi- lidade que não tem contraparte na forma de Lagrange. Isso também simplifica a aplicação do Princípio do Erro (utilizado no Exemplo 3 para Lagrange). Os detalhes dessas idéias são os seguintes: Consideremos que pn–1(x) seja o (n – 1)-ésimo polinômio de Newton (cuja forma determinaremos); assim, pn-1(x0) = f0, pn-1(x1) = f1,• • •, pn-1(xn-1) = fn–1. Além disso, escrevamos o n-ésimo polinômio de Newton como (6) pn(x) pn1(x) gn(x); donde (6 ) gn(x) pn(x) pn1(x), Aqui, gn(x) deve ser determinado de modo que pn(x0) ƒ0, pn(x1) ƒ1, • • • , pn(xn) ƒn. Como pn e pn-1 concordam em x0,• • •, xn-1, vemos que gn é zero aí. Além disso, em geral gn será um polinômio de grau, no máximo, igual a n – 1. Logo, gn deve ter a forma (6 ) gn(x) an(x x0)(x x1) • • • (x xn1). Determinemos a constante an. Para isto, façamos x = xn e resolvamos algebricamente (6 ) para an. Substituindo gn(xn) de acordo com (6 ) e usando pn(xn) = fn, vemos que isso resulta em (7) an ƒn pn1(xn) (xn x0)(xn x1) • • • (xn xn1) . Escrevamos ak em vez de an e mostremos que ak é igual à k-ésima diferença dividida, que é recursivamente denotada e definida do seguinte modo: a1 ƒ[x0, x1] ƒ1 ƒ0 x1 x0 a2 ƒ[x0, x1, x2] ƒ[x1, x2] ƒ[x0, x1] x2 x0 e, em geral, (8) ak ƒ[x0, • • • , xk] ƒ[x1, • • • , xk] ƒ[x0, • • • , xk1] xk x0 . PROVA Se n = 1, então pn1(xn) p0(x1) ƒ0, pois p0(x) é constante e igual a f0, o valor de f(x) em x0. Logo, (7) fornece a1 ƒ1 p0(x1) x1 x0 ƒ1 ƒ0 x1 x0 ƒ[x0, x1], com (6) e (6 ) fornecendo o polinômio de interpolação de Newton do primeiro grau p1(x) ƒ0 (x x0)ƒ[x0, x1]. Se n = 2, então este p1 e (7) dão a2 ƒ2 p1(x2) (x2 x0)(x2 x1) ƒ2 ƒ0 (x2 x0)ƒ[x0, x1] (x2 x0)(x2 x1) ƒ[x0, x1, x2] 24 Parte E • Análise Numérica Em (10), finalmente fazemos x = x0 + rh. Então, x x0 rh, x x1 (r 1)h, visto que x1 – x0 = h, e assim por diante. Com isto e com (13), a fórmula (10) passa a ser a fórmula de interpolação da diferença progressiva de Newton (ou de Gregory2–Newton) ƒ(x) pn(x) n s0 rs sƒ0 (x x0 rh, r (x x0) /h)(14) ƒ0 rƒ0 r(r 1) 2! 2ƒ0 • • • r(r 1) • • • (r n 1) n! nƒ0 onde os coeficientes binomiais na primeira linha são definidos por (15) r0 1, r s r(r 1)(r 2) • • • (r s 1) s! (s 0, inteiro) e s! 1 • 2 • • • s. Erro. De (5), obtemos, com x x0 rh, x x1 (r 1)h etc., (16) en(x) ƒ(x) pn(x) hn 1 (n 1)! r(r 1) • • • (r n)ƒ (n 1)(t) com t conforme caracterizado em (5). A fórmula (16) é uma fórmula exata para o erro, porém envolve a incógnita t. No Exemplo 5 (a seguir), mos- traremos como usar (16) para obter uma estimativa de erro e um intervalo dentro do qual deve se situar o valor verdadeiro de f(x). Comentários sobre a Precisão. (A) A ordem de magnitude do erro en(x) é aproximadamente igual à da diferença seguinte não utilizada em pn(x). (B) Devemos escolher x0,• • •, xn tais que o valor de x no qual a interpolação é feita esteja tão bem centrado entre x0,• • •, xn quanto possível. A razão para (A) é que, em (16), ƒn 1(t) n 1ƒ(t) hn 1 , r(r 1) • • • (r n) 1 2 • • • (n 1) 1 se r 1 (e, na verdade, para qualquer r, contanto que não extrapolemos). A razão para (B) é que r(r 1) • • • (r n) torna-se o menor valor para essa escolha. EXEMPLO 5 Fórmula da Diferença Progressiva de Newton. Estimativa do Erro Calcule cosh 0,56 a partir de (14) e dos quatro valores na tabela a seguir, e faça a estimativa do erro. j xj ƒj cosh xj ƒj 2ƒj 3ƒj 0 0,5 1,127 626 0,057 839 1 0,6 1,185 465 0,011 865 0,069 704 0,000 697 2 0,7 1,255 169 0,012 562 0,082 266 3 0,8 1,337 435 Solução. Calculemos as diferenças progressivas conforme mostra a tabela. Os valores de que precisamos estão circulados. Em (14), temos r (0,56 0,50)/0,1 0,6, de modo que (14) fornece 2 JAMES GREGORY (1638–1675), matemático escocês e catedrático nas universidades de St. Andrews e Edimburgo. Os símbolos em (14) e 2 nada têm a ver com o laplaciano. Capítulo 19: Métodos Numéricos em Geral 25 cosh 0,56 1,127 626 0,6 0,057 839 0,6(0,4) 2 0,011 865 0,6(0,4)(1,4) 6 0,000 697 1,127 626 0,034 703 0,001 424 0,000 039 1,160 944. Estimativa do erro. De (16), visto que a quarta derivada é cosh(4) t = cosh t, e3(0,56) 0,14 4! 0,6(0,4)(1,4)(2,4) cosh t A cosh t, onde A 0,000 003 36 e 0,5 t 0,8. Não conhecemos t, mas obtemos uma desigualdade tomando os valores máximo e mínimo de cosh t nesse intervalo: A cosh 0,8 e3(0,62) A cosh 0,5. Como ƒ(x) p3(x) e3(x), isto fornece p3(0,56) A cosh 0,8 cosh 0,56 p3(0,56) A cosh 0,5. Os valores numéricos são 1,160 939 cosh 0,56 1,160 941. O valor exato em 6D é cosh 0,56 = 1,160 941, que está situado dentro desses limites. Nem sempre tais limites são tão estreitos. Além disso, não consideramos os erros de arredondamento, que dependerão do número de operações. Este exemplo também explica o nome “fórmula da diferença progressiva”: vemos que as diferenças na fórmula inclinam-se para a frente na tabela de diferenças. Espaçamento Igual: Fórmula da Diferença Regressiva de Newton Ao invés das diferenças de inclinação progressiva, podemos também empregar as diferenças de inclinação regres- siva. A tabela de diferenças continua a mesma de antes (os mesmos números, na mesma posição), exceto por uma alteração muito inofensiva do número subscrito em questão j (que explicaremos no Exemplo 6, a seguir). Não obstante, puramente por razões de conveniência, é padrão apresentar um segundo nome e uma segunda notação para as diferenças, conforme se segue. Definimos a primeira diferença regressiva de f em xj por ƒj ƒj ƒj1, a segunda diferença regressiva de f em xj por 2ƒj ƒj ƒj1, e, continuando dessa maneira, a k-ésima diferença regressiva de f em xj por (17) kƒj k1ƒj k1ƒj1 (k 1, 2, • • •). Uma fórmula semelhante a (14) mas que envolve diferenças regressivas é a fórmula de interpolação de diferença regressiva de Newton (ou de Gregory–Newton). (18) ƒ(x) pn(x) n s0 r s 1s sƒ0 (x x0 rh, r (x x0) /h) ƒ0 rƒ0 r(r 1) 2! 2ƒ0 • • • r(r 1) • • • (r n 1) n! nƒ0. EXEMPLO 6 Interpolações Progressivas e Regressivas de Newton Calcule um valor em 7D da função de Bessel J0(x) para x = 1,72 a partir dos valores da tabela a seguir, usando: (a) a fórmula progressiva de Newton (14) e (b) a fórmula regressiva de Newton (18). 26 Parte E • Análise Numérica jprogr jregr xj J0(xj) 1a Deriv. 2a Deriv. 3a Deriv. 0 3 1,7 0,397 9849 0,057 9985 1 2 1,8 0,339 9864 0,000 1693 0,058 1678 0,000 4093 2 1 1,9 0,281 8186 0,000 2400 0,057 9278 3 0 2,0 0,223 8908 Solução. O cálculo das diferenças é o mesmo em ambos os casos. Somente a sua notação é que difere. (a) Progressiva. Em (14), temos r (1,72 1,70)/0,1 0,2 e j vai de 0 a 3 (veja a primeira coluna). Em cada coluna, precisamos do primeiro número dado e, assim, (14) fornece J0(1,72) 0,397 9849 0,2(0,057 9985) 0,2(0,8) 2 (0,000 1693) 0,2(0,8)(1,8) 6 0,000 4093 0,397 9849 0,011 5997 0,000 0135 0,000 0196 0,386 4183, que é um valor exato em 6D, com o valor exato em 7D sendo 0,386 4185. (b) Regressiva. Para (18), usamos o j mostrado na segunda coluna e, em cada coluna, o último número. Como r (1,72 2,00)/0,1 2,8, obtemos então de (18) J0(1,72) 0,223 8908 2,8(0,057 9278) 2,8(1,8) 2 0,000 2400 2,8(1,8)(0,8) 6 0,000 4093 0,223 8908 0,162 1978 0,000 6048 0,000 2750 0,386 4184. Notação da Diferença Central Esta é uma terceira notação para as diferenças. A primeira diferença central de f(x) em xj é definida por dƒj ƒj 1/2 ƒj1/2 e a k-ésima diferença central de f(x) em xj, por (19) dkƒj dk1ƒj 1/2 dk1ƒj1/2 ( j 2, 3, • • •). Assim, nesta notação, uma tabela de diferenças, por exemplo, para f-1, f0, f1, f2, tem a seguinte aparência: x1 x0 x1 x2 ƒ1 ƒ0 ƒ1 ƒ2 dƒ1/2 dƒ1/2 dƒ3/2 d2ƒ0 d2ƒ1 d3ƒ1/2 Utilizam-se as diferenças centrais na derivação numérica (Seção 19.5), em equações diferenciais (Capítulo 21) e em fórmulas centralizadas de interpolação (p. ex., a fórmula de Everett no Projeto de Equipe 22). Estas são fórmulas que utilizam valores de funções localizados “simetricamente” em ambos os lados do ponto de interpo- lação x. Esses valores estão disponíveis próximo do meio de uma tabela dada, em que as fórmulas centralizadas de interpolação tendem a fornecer resultados melhores que os das fórmulas de Newton, as quais não possuem essa propriedade de “simetria”. P R O B L E M A S P R O P O S T O S 1 9 . 3 1. (Interpolação linear) Calcule p1(x) no Exemplo 1. A partir dele, calcule ln 9,4 ≈ p1(9,4). 2. Faça a estimativa do erro no Problema 1, utilizando (5). 3. (Interpolação quadrática) Calcule o polinômio de Lagrange p2(x) para os valores em 4D da função gama [(24), Apêndice 3.1] (1,00) 1,0000, (1,02) 0,9888, (1,04) 0,9784 e, a partir Capítulo 19: Métodos Numéricos em Geral 29 Os qj’s mais simples seriam polinômios lineares. Entretanto, a curva de uma função linear contínua por inter- valos apresenta arestas, e teria pouco interesse geral — imagine, por exemplo, a modelagem da estrutura de um automóvel ou de um navio. Consideraremos os splines cúbicos por serem os mais importantes nas aplicações. Por definição, um spline cúbico g(x) interpolando dados fornecidos (x0, f0), • • • , (xn, fn) é uma função contínua no intervalo a = x0 x xn = b, que tem a primeira e a segunda derivativas contínuas e que satisfaz à condição de interpolação (2); além disso, entre nós adjacentes, g(x) é dada por um polinômio qj(x) de grau menor ou igual a 3. Afirmamos que esse spline cúbico existe. E se, além de (2), também exigirmos que (3) g (x0) k0, g (xn) kn (direções tangentes a g(x) dadas nas duas extremidades do intervalo a x b), temos então um spline cúbico unicamente determinado. Este é o conteúdo do seguinte teorema da existência e unicidade, cuja prova também sugerirá a determinação real dos splines. (A Condição (3) será discutida após a demonstração.) TEOREMA 1 Existência e Unicidade de Splines Cúbicos Consideremos os pontos (x0, f0), (x1, f1), • • • , (xn, fn) com os xj arbitrariamente espaçados [veja (1)] e os valores dados fj = f(xj), j = 0, 1, • • • , n. Consideremos também que k0 e kn sejam quaisquer números dados. Então, existe um e somente um spline cúbico g(x) que corresponde a (1) e satisfaz a (2) e (3). PROVA Por definição, em cada subintervalo Ij dado por xj x xj+1, o spline g(x) precisa concordar com um polinômio qj(x) de grau não superior a 3, de modo que (4) qj(xj) ƒ(xj), qj(xj 1) ƒ(xj 1) ( j 0, 1, • • • , n 1). Para as derivadas, escrevemos; (5) q j(xj) kj, q j(xj 1) kj 1 ( j 0, 1, • • • , n 1) com k0 e kn dados e k1, • • • , kn-1 a serem determinados posteriormente. As Equações (4) e (5) são quatro condições a serem satisfeitas por cada qj(x). Por cálculo direto e empregando a notação (6*) cj 1 hj 1 xj+1 – xj ( j 0, 1, • • • , n 1) podemos verificar que o único polinômio cúbico qj(x) (j = 0, 1, • • • , n –1) que satisfaz a (4) e (5) é (6) ƒ(xj)cj 2(x xj 1) 2[1 2cj(x xj)] ƒ(xj 1)cj 2(x xj) 2[1 2cj(x xj 1)] kjcj 2(x xj)(x xj 1) 2 kj 1cj 2(x xj) 2(x xj 1). qj(x) Derivando duas vezes, obtemos (7) q j(xj) 6cj2ƒ(xj) 6cj2ƒ(xj 1) 4cjkj 2cjkj 1 (8) q j(xj 1) 6cj2ƒ(xj) 6cj2ƒ(xj 1) 2cjkj 4cjkj 1. Por definição, g(x) possui derivadas segundas contínuas. Isto fornece as condições q j1(xj) q j(xj) ( j 1, • • • , n 1). Se usarmos (7) e também (8) com a substituição de j por j – 1, essas n – 1 equações tornam-se (9) cj1kj1 2(cj1 cj)kj cjkj 1 3[c2j1ƒj cj2 ƒj 1] onde ƒj ƒ(xj) ƒ(xj1) e ƒj 1 ƒ(xj 1) ƒ(xj) e j 1, • • • , n 1, como antes. Esse sistema linear de n – 1 equações tem uma solução única k1, • • • , kn-1, visto que a matriz dos coeficientes é estritamente dominante diagonalmente [ou seja, em cada linha, o elemento diagonal (positivo) é maior do que a soma dos outros elementos 30 Parte E • Análise Numérica (positivos)]. Logo, o determinante da matriz não pode ser zero (como decorre do Teorema 3 na Seção 20.7), de modo que podemos determinar valores únicos k1, • • • , kn-1 da primeira derivada de g(x) nos nós. Isto prova o teorema. As Demandas de Armazenamento e de Tempo na resolução de (9) são modestas, já que a matriz de (9) é esparsa (tem poucos elementos não-nulos) e tridiagonal (pode ter elementos não-nulos somente na diagonal e nas duas “paralelas” adjacentes, acima e abaixo dela). A ocorrência de dominância dispensa a necessidade de pivotação (Seção 7.3). Isto torna os splines eficientes na resolução de problemas grandes, com milhares de nós ou mais. Sobre referências em literatura e alguns comentários críticos, veja a American Mathematical Monthly 105 (1998), 929–941. A condição (3) inclui as condições fixas (10) g (x0) ƒ (x0), g (xn) ƒ (xn), nas quais são dadas as direções f (x0) e f (xn) tangentes às extremidades. Outras condições de interesse prático são as condições livres ou naturais. (11) g (x0) 0, g (xn) 0 (geometricamente: curvatura zero nas extremidades, como ocorre na régua do projetista), que fornecem um spline natural. Estes nomes são justificados pela Fig. 290, em Problemas Propostos 12.3. Determinação de Splines. Considere que k0 e kn sejam dados. Obtenha k1, • • • , kn-1 resolvendo o sistema linear (9). Lembre-se de que o spline g(x) a ser encontrado consiste em n polinômios cúbicos q0, • • • , qn-1. Escrevemos esses polinômios na forma (12) qj(x) aj0 aj1(x xj) aj2(x xj)2 aj3(x xj)3 onde j = 0, • • • , n – 1. Utilizando a fórmula de Taylor, obtemos aj0 qj(xj) ƒj por (2), aj1 q j(xj) kj por (5), (13) aj2 1 2 q j(xj) 3 h2 (ƒj 1 ƒj) 1 h (kj 1 2kj) por (7), aj3 1 6 qj(xj) 2 h3 (ƒj ƒj 1) 1 h2 (kj 1 kj) com aj3 sendo obtido calculando-se q j(xj 1) a partir de (12) e igualando o resultado a (8), ou seja: q j(xj 1) 2aj2 6aj3h 6 h2 (ƒj ƒj 1) 2 h (kj 2kj 1), e então subtraindo disto 2aj2, conforme dado em (13), e simplificando. Note que, para os nós eqüidistantes de distância hj = h podemos escrever cj = c = 1/h em (6*) e obter, a partir de (9), simplesmente (14) kj1 4kj kj 1 3 h (ƒj 1 ƒj1) ( j 1, • • • , n 1). EXEMPLO 1 Interpolação Spline. Nós Eqüidistantes Interpole f(x) = x4 no intervalo –1 x 1 usando o spline cúbico g(x) correspondente aos nós x0 = –1, x1 = 0, x2 = 1 e satisfazendo às condições fixas g (–1) = f (–1), g (1) = f (1). Solução. Em nossa notação-padrão, os dados fornecidos são ƒ0 ƒ(1) 1, ƒ1 ƒ(0) 0, ƒ2 ƒ(1) 1. Temos que h = 1 e n = 2, de modo que nosso spline consiste em n = 2 polinômios q0(x) a00 a01(x 1) a02(x 1) 2 a03(x 1) 3 (1 x 0), q1(x) a10 a11x a12x 2 a13x 3 (0 x 1). Capítulo 19: Métodos Numéricos em Geral 31 Determinamos o valor de kj a partir de (14) (eqüidistância!) e então os coeficientes do spline a partir de (13). Como n = 2, o sistema (14) é uma equação única (com j = 1 e h = 1) k0 4k1 k2 3(ƒ2 ƒ0). Aqui, f0 = f2 = 1 (o valor de x4 nas extremidades) e k0 = –4, k2 = 4, com os valores da derivada 4x3 nas extremidades sendo –1 e 1. Logo, 4 4k1 4 3(1 1) 0, k1 0. De (13), podemos agora obter os coeficientes de q0, ou seja, a00 ƒ0 1, a01 k0 4, e 5. 2. 3(0 1) (0 8) 2(1 0) (0 4) 1 3 2 (ƒ1 ƒ0) 1 1 (k1 2k0) 1 2 3 (ƒ0 ƒ1) 1 1 2 (k1 k0) a02 a03 De modo semelhante, para os coeficientes de q1, obtemos, a partir de (13), os valores a10 ƒ1 0, a11 k1 0 e 3(1 0) (4 0) 1 2(0 1) (4 0) 2. 2k1) k1) 3(ƒ2 ƒ1) (k2 2(ƒ1 ƒ2) (k2 a12 a13 Isso fornece os polinômios nos quais o spline g(x) consiste, ou seja, g(x) 1 x 0.0 x 1. se se q0(x) 1 4(x 1) 5(x 1) 2 2(x 1)3 x2 2x3 q1(x) x 2 2x3 A Fig. 433 mostra f(x) e este spline. Você percebe que poderíamos ter poupado quase a metade do nosso trabalho se tivéssemos usado a simetria do problema? 1 –1 1 f(x) g(x) x Fig. 433. Função f(x) = x4 e spline cúbico g(x) do Exemplo 1 EXEMPLO 2 Spline Natural. Nós Arbitrariamente Espaçados Encontre uma aproximação por spline e uma aproximação polinomial para a curva de secção transversal do Santuário do Livro em Jerusalém, que tem forma circular e está mostrado na Fig. 434. –6 –5 –4 –3 –2 –1 0 1 1 2 3 Fig. 434. Santuário do Livro em Jerusalém (Arquitetos F. Kissler e A. M. Bartus) Solução. Treze pontos, distribuídos de forma aproximadamente igual ao longo do contorno (e não ao longo do eixo x!), fornecem os seguintes números: 34 Parte E • Análise Numérica Regra do Retângulo. Regra do Trapézio Os métodos de integração numérica são obtidos conseguindo-se aproximações do integrando f por meio de funções que podem ser integradas facilmente. A fórmula mais simples, correspondente à regra do retângulo, é obtida se subdividirmos o intervalo de integração a x b em n intervalos de igual comprimento h = (b – a)/n e, em cada subintervalo, obtivermos uma aproximação de f pela constante f(xj*), que é o valor de f no ponto médio xj* do j-ésimo subintervalo (Fig. 438). Então, obtemos uma aproximação de f por meio de função degrau (função constante por intervalos), os n retângulos na Fig. 438 têm as áreas f(x1*)h, • • •, f(xn*)h, e a regra do retângulo é (1) J b a ƒ(x) dx h[ƒ(x1*) ƒ(x2*) • • • ƒ(xn*)] h b a n . A regra do trapézio é geralmente mais precisa. Obtêmo-la se tomarmos a mesma subdivisão de antes e apro- ximarmos f por uma linha quebrada de segmentos (cordas) com as extremidades [a, f(a)], [x1, f(x1)], • • •, [b, f(b)] na curva de f (Fig. 439). Então a área abaixo da curva de f entre a e b é aproximada por n trapézios de áreas 1 2 [ƒ(a) ƒ(x1)]h, 12[ƒ(x1) ƒ(x2)]h, • • • , 1 2 [ƒ(xn1) ƒ(b)]h. y xa b y = f(x) R y xa b y = f(x) x 1 * x 2 * x n *Fig. 437. Interpretação geométrica de uma integral definida Fig. 438. Regra do retângulo y xa b y = f(x) x 1 x 2 x n – 1 Fig. 439. Regra do trapézio Fazendo a soma desses trapézios, obtemos a regra do trapézio (2) J b a ƒ(x) dx h[1_2ƒ(a) ƒ(x1) ƒ(x2) • • • ƒ(xn1) 1_2ƒ(b)] onde h = (b – a)/n, como em (1). Os xj’s e a e b são chamados de nós. EXEMPLO 1 Regra do Trapézio Avalie J 1 0 ex 2 dx usando (2) com n = 10. Solução. J ≈ 0,1(0,5 1,367 879 + 6,778 167) = 0,746 211 da Tabela 19.3. Tabela 19.3 Cálculos do Exemplo 1 j xj xj2 e xj 2 0 0 0 1,000 000 1 0,1 0,01 0,990 050 2 0,2 0,04 0,960 789 3 0,3 0,09 0,913 931 4 0,4 0,16 0,852 144 5 0,5 0,25 0,778 801 6 0,6 0,36 0,697 676 7 0,7 0,49 0,612 626 8 0,8 0,64 0,527 292 9 0,9 0,81 0,444 858 10 1,0 1,00 0,367 879 Somas 1,367 879 6,778 167 Capítulo 19: Métodos Numéricos em Geral 35 Limites e Estimativa de Erro para a Regra do Trapézio Uma estimativa de erro para a regra do trapézio pode ser obtida de (5) na Seção 19.3 com n = 1, usando-se integração como se segue. Para um único subintervalo, temos ƒ(x) p1(x) (x x0)(x x1) ƒ ( t) 2 com um valor apropriado de t dependente de x, situado entre x0 e x1. A integração em x de a = x0 até x1 = x0 + h fornece x0 x0 h ƒ(x) dx h 2 [ƒ(x0) ƒ(x1)] x0 x0 h (x x0)(x x0 h) ƒ (t(x)) 2 dx. Fazendo x – x0 = v e aplicando o teorema do valor médio do cálculo integral, que podemos utilizar, visto que (x – x0)(x – x0 – h) não muda de sinal, encontramos que o lado direito se iguala a (3*) h 0 v(v h) dv ƒ ( t) 2 h 3 3 h3 2 ƒ ( t) 2 h3 12 ƒ ( t) onde t é um valor (adequado e desconhecido) situado entre x0 e x1. Este é o erro para a regra do trapézio com n = 1, freqüentemente chamado de erro local. Logo, o erro e de (2) com n qualquer é o somatório dessas contribuições provindas dos n subintervalos; como h = (b – a)/n, nh3 = n(b – a)3/n3, e (b – a)2 = n2h2, obtemos (3) e (b a)3 12n2 ƒ (t̂) (b a) 12 h2ƒ (t̂) com t̂ (adequado e desconhecido) situado entre a e b. Devido à Equação (3), a regra do trapézio é também escrita como (2*) J b a ƒ(x) dx h [12ƒ(a) ƒ(x1) • • • ƒ(xn1) 12ƒ(b)] (b a) 12 h2ƒ (t̂). Os Limites de Erro são agora obtidos tomando-se o maior valor para f , digamos, M2, e o menor valor, M2*, no intervalo de integração. Então, (3) fornece (note que K é negativo) (4) KM2 e KM2* onde K (b a)3 12n2 b a 12 h2. A Estimativa do Erro Reduzindo-se h à Metade é aconselhável se h for muito complicado ou desconhecido, por exemplo, no caso de os dados serem experimentais. Então, podemos aplicar o Princípio do Erro da Seção 19.1. Ou seja, fazemos o cálculo usando (2), primeiramente com h, obtendo, digamos, J = Jh + eh, e então com 1 2 h, obtendo J = Jh/2 + eh/2. Ora, se em (3) substituirmos h2 por (12h) 2, o erro fica multiplicado por 1 4 . Logo, eh/2 ≈ 14eh (não exatamente igual porque t̂ pode diferir). Juntos, Jh/2 + eh/2 = Jh + eh ≈ Jh + 4eh/2. Portanto, Jh/2 – Jh = (4 – 1)eh/2. Dividindo por 3, obtemos a fórmula do erro para Jh/2 (5) eh/2 1 3 (Jh/2 Jh). EXEMPLO 2 Estimativa do Erro para a Regra do Trapézio Usando-se (4) e (5) Estime o erro do valor aproximado no Exemplo 1 utilizando (4) e (5). Solução. (A) Limites do erro por (4). Por derivação, f (x) = 2(2x2 – 1)ex 2 . Além disso, f (x) > 0 se 0 < x < 1, de modo que os valores mínimo e máximo ocorrem nas extremidades do intervalo. Calculamos M2 = f (1) = 0,735 759 e M2* = f (0) = –2. Além disso, K = –1/1200, e (4) fornece –0,000 614 e 0,001 667. Logo, o valor exato de J deve se situar entre 0,746 211 – 0,000 614 = 0,745 597 e 0,746 211 + 0,001 667 = 0,747 878. Na verdade, J = 0,746 824, com uma exatidão de 6D. 36 Parte E • Análise Numérica (B) Estimativa do erro por (5). Jh = 0,746 211 no Exemplo 1. Além disso, Jh/2 0,05 19 j 1 e ( j/20) 2 1 2 (1 0,367879) 0,746671. Logo, eh/2 = 1 3(Jh/2 – Jh) = 0,000 153 e Jh/2 + eh/2 = 0,746 824, com uma exatidão de 6D. Regra de Integração de Simpson Se uma aproximação constante por intervalos de f nos levou à regra do retângulo (1) e uma aproximação linear por intervalos nos levou à regra do trapézio (2), agora uma aproximação quadrática por intervalos nos levará à regra de Simpson, que é de grande importância prática devido ao fato de ser suficientemente exata para a maioria dos problemas, embora ainda suficientemente simples. Para obtermos a regra de Simpson, dividimos o intervalo de integração a x b em um número par de subintervalos iguais, digamos, em n = 2m subintervalos de comprimento h = (b – a)/2m, com as extremidades x0 (= a), x1, • • • , x2m-1, x2m (= b); veja a Fig. 440. Agora tomamos os dois primeiros subintervalos e obtemos uma aproximação de f(x) no intervalo x0 x x2 = x0 + 2h usando o polinômio de Lagrange p2(x) que passa pelos pontos (x0, f0), (x1, f1), (x2, f2), onde fj = f(xj). De (3) da Seção 19.3, obtemos (6) p2(x) (x x1)(x x2) (x0 x1)(x0 x2) ƒ0 (x x0)(x x2) (x1 x0)(x1 x2) ƒ1 (x x0)(x x1) (x2 x0)(x2 x1) ƒ2. Os denominadores em (6) são 2h2, –h2 e 2h2, respectivamente. Fazendo s = (x – x1)/h, temos x x1 sh, x x0 x (x1 h) (s 1)h x x2 x (x1 h) (s 1)h e obtemos p2(x) 12s (s 1)ƒ0 (s 1)(s 1)ƒ1 1 2 (s 1)sƒ2. Agora integramos em x de x0 a x2. Isto corresponde a integrar em s de –1 a 1. Como dx = h ds, o resultado é (7*) x0 x2 ƒ(x) dx x0 x2 p2(x) dx h 1 3 ƒ0 4 3 ƒ1 1 3 ƒ2 . Uma fórmula similar se verifica para os dois próximos subintervalos de x2 a x4, e assim por diante. Somando todas essas m fórmulas, obtemos a regra de Simpson4 (7) b a ƒ(x) dx h 3 (ƒ0 4ƒ1 2ƒ2 4ƒ3 • • • 2ƒ2m2 4ƒ2m1 ƒ2m), onde h = (b – a)/(2m) e fj = f(xj). A Tabela 19.4 apresenta um algoritmo para a regra de Simpson. y xa bx1 x2 x3 x4 x2m–2 x2m–1 Primeira parábola Última parábola Segunda parábola y = f(x) Fig. 440. Regra de Simpson 4THOMAS SIMPSON (1710–1761), matemático autodidata inglês, autor de diversos livros-texto populares. A regra de Simpson foi usada muito mais cedo por Torricelli, Gregory (em 1668) e Newton (em 1676). Capítulo 19: Métodos Numéricos em Geral 39 (10) eh/2 1 15 (Jh/2 Jh). Jh é obtido usando-se h; Jh/2 é obtido usando-se 12h, e eh/2 é o erro de Jh/2. Obtenção. Em (5), tínhamos 1 3 como o recíproco de 3 = 4 – 1 e 1 4 = (1 2 )2 resultante de h2 em (3) substituindo-se h por 1 2 h. Em (10), temos 1 15 como o recíproco de 15 = 16 – 1 e 1 16 = (1 2 )4 resultante de h4 em (8) substituindo-se h por 1 2 h. EXEMPLO 5 Estimativa de Erro para a Regra de Simpson pela Redução de h à Metade Integre ƒ(x) 14px 4 cos 14px de 0 a 2 com h = 1 e aplique (10). Solução. O valor exato 5D da integral é J = 1,25953. A regra de Simpson fornece Jh 1 3[ƒ(0) 4ƒ(1) ƒ(2)] 1 3(0 4 0,555360 0) 0,740480, Jh/2 1 6 ƒ(0) 4ƒ 12 2ƒ(1) 4ƒ 3 2 ƒ(2) 1 6 [0 4 0,045351 2 0,555361 4 1,521579 0] 1,22974. Logo (10) fornece eh/2 = 1 15(1,22974 – 0,74048) = 0,032617 e, assim, J ≈ Jh/2 + eh/2 = 1,26236, com um erro de –0,00283, que é menor em módulo do que 110 do erro 0,02979 de Jh/2. Portanto, o uso de (10) foi bastante proveitoso. Integração Adaptativa A idéia é adaptar o passo h à variabilidade de f(x). Ou seja, onde f varia apenas um pouco, podemos proceder em grandes passos sem provocar um erro substancial na integral, mas onde f varia rapidamente, temos de efetuar pequenos passos a fim de, em todos os lugares, permanecermos perto o bastante da curva de f. A mudança de h é feita sistematicamente, usualmente dividindo-se h pela metade, e automaticamente (não “de forma manual”) dependendo do tamanho do erro (estimado) ao longo de um subintervalo. O subintervalo é dividido pela metade se o erro correspondente for ainda muito grande, ou seja, maior do que uma dada tolerância TOL (máximo erro absoluto admissível), ou não é dividido pela metade se o erro for menor ou igual à TOL. A adaptação é uma das técnicas típicas utilizadas nos modernos softwares. Em conexão com a integração, ela pode ser aplicada a vários métodos. Explicá-la-emos aqui no contexto da regra de Simpson. Na Tabela 19.6, um asterisco indica que se atingiu TOL para o subintervalo em questão. EXEMPLO 6 Integração Adaptativa com a Regra de Simpson Integre ƒ(x) 14px 4 cos 14px de x = 0 a x = 2, usando a integração adaptativa e a regra de Simpson, com TOL[0,2] = 0,0002. Solução. A Tabela 19.6 mostra os cálculos. A Fig. 441 mostra o integrando f(x) e os intervalos adaptados utilizados. Os dois primeiros intervalos ([0, 0,5], [0,5, 1,0]) têm um comprimento de 0,5, logo, h = 0,25 [pois usamos 2m = 2 subintervalos na regra de Simpson (7**)]. Os dois intervalos seguintes ([1,00, 1,25], [1,25, 1,50]) têm um comprimento de 0,25 (logo, h = 0,125) e os quatro últimos intervalos têm um comprimento de 0,125. Cálculos amostrais. Para 0,740480, veja o Exemplo 5. A fórmula (10) fornece (0,123716 – 0,122794)/15 = 0,000061. Note que 0,123716 se refere a [0, 0,5] e [0,5, 1], de modo que precisamos subtrair o valor correspondente a [0,1] na linha anterior etc. A TOL[0,2] = 0,0002 fornece 0,0001 para os subintervalos de comprimento 1, fornece 0,00005 para os subintervalos de comprimento 0,5 etc. O valor da integral que se obtém é a soma dos valores marcados por um asterisco (para os quais a estimativa de erro torna-se menor do que a TOL). Isso nos dá J 0,123716 0,528895 0,388263 0,218483 1,25936. O valor exato em 5D é J = 1,25953. Logo, o erro é 0,00017. Isto é cerca de 1/200 do valor absoluto do erro no Exemplo 5. O cálculo mais longo que agora fizemos produziu um resultado muito melhor. f(x) x0 1,5 1,0 0,5 0,50 1 1,5 2 Fig. 441. Integração adaptativa do Exemplo 6 40 Parte E • Análise Numérica Tabela 19.6 Cálculos do Exemplo 6 Intervalo Integral Erro (10) TOL Comentário [0, 2] 0,740480 0,0002 [0, 1] 0,122794 [1, 2] 1,10695 Soma 1,22974 0,032617 0,0002 Dividir novamente [0,0, 0,5] 0,004782 [0,5, 1,0] 0,118934 Soma 0,123716* 0,000061 0,0001 TOL alcançada [1,0, 1,5] 0,528176 [1,5, 2,0] 0,605821 Soma 1,13300 0,001803 0,0001 Dividir novamente [1,00, 1,25] 0,200544 [1,25, 1,50] 0,328351 Soma 0,528895* 0,000048 0,00005 TOL alcançada [1,50, 1,75] 0,388235 [1,75, 2,00] 0,218457 Soma 0,606692 0,000058 0,00005 Dividir novamente [1,500, 1,625] 0,196244 [1,625, 1,750] 0,192019 Soma 0,388263* 0,000002 0,000025 TOL alcançada [1,750, 1,875] 0,153405 [1,875, 2,000] 0,065078 Soma 0,218483* 0,000002 0,000025 TOL alcançada Fórmulas de Integração de Gauss. Máximo Grau de Precisão As fórmulas de integração que discutimos até agora utilizam valores da função em valores de x (nós) predetermi- nados (eqüidistantes) e fornecem resultados exatos para polinômios que não excedem um certo grau [chamado de grau de precisão; veja depois (9)]. Porém, podemos obter fórmulas de integração muito mais exatas do seguinte modo. Fazemos (11) 1 1 ƒ(t) dt n j1 Ajƒj [ƒj ƒ(tj)] com n fixo, e t = 1 obtido de x = a, b fazendo x = 1 2 [a(t – 1) + b(t + 1)]. Então, determinamos os n coeficientes A1, • • • , An e n nós t1, • • • , tn, de modo que (11) forneça resultados exatos para polinômios de um grau k tão alto quanto possível. Como n + n = 2n é o número de coeficientes de um polinômio de grau 2n – 1, segue-se que k 2n – 1. Gauss mostrou que é possível atingir uma precisão para polinômios de grau não excedendo 2n – 1 (ao invés de n – 1, como ocorria com os nós predeterminados), e ele forneceu a localização de tj (= o j-ésimo zero do polinômio de Legendre Pn na Seção 5.3) e dos coeficientes Aj que dependem de n, mas não de f(t), e que são obtidos usando-se o polinômio de interpolação de Lagrange, conforme mostra a Ref. [E5] no Apêndice 1. Com esses tj e Aj, a fórmula (11) é chamada de fórmula de integração de Gauss ou fórmula da quadratura de Gauss. Seu grau de precisão é 2n – 1, como dizemos há pouco. A Tabela 19.7 fornece os valores necessários para n = 2, • • • , 5. (Para n maior, veja a Ref. [RG1] no Apêndice 1.) EXEMPLO 7 Fórmula de Integração de Gauss com n = 3 Avalie a integral no Exemplo 3 pela fórmula de integração de Gauss (11) com n = 3. Solução. Temos que converter nossa integral de 0 a 1 em uma integral de –1 a 1. Fazemos então x = 12 (t + 1). Logo, dx = 1 2 dt, e (11) com n = 3 e os valores acima dos nós e os coeficientes fornecem Capítulo 19: Métodos Numéricos em Geral 41 1 0 exp ( x2) dx 1 1 exp ( (t 1)2) dt exp 1 ) 2 exp exp 1 2 0,746 815351459148935145912 1 4 1 2 (valor exato em 6D: 0,746 825), que é quase tão exato quanto o resultado obtido pela regra de Simpson no Exemplo 3, onde era necessário um número muito maior de operações aritméticas. Usando 3 valores da função (como neste exemplo) e a regra de Simpson, obteríamos 16 (1 + 4e–0,25 + e–1) = 0,747 180, com um erro mais de 30 vezes superior ao da integração de Gauss. EXEMPLO 8 Fórmula de Integração de Gauss com n = 4 e 5 Integre ƒ(x) 14px 4 cos 14px de 0 a 2 por Gauss. Compare com a integração adaptativa no Exemplo 6 e comente. Solução. x = t + 1 fornece ƒ(t) 14p(t 1)4 cos ( 1 4p(t 1)), conforme necessário em (11). Para n = 4, calculamos (6S) J A1ƒ1 • • • A4ƒ4 A1(ƒ1 ƒ4) A2(ƒ2 ƒ3) 0,347855(0,000290309 1,02570) 0,652145(0,129464 1,25459) 1,25950. O erro é 0,00003, pois J = 1,25953 (6S). Calculando com 10S e n = 4, obtemos o mesmo resultado; logo, o erro se deve à fórmula, não ao arredondamento. Para n = 5 e 10S, temos J ≈ 1,25952 6185, um valor maior pela quantidade 0,00000 0250, pois J = 1,25952 5935 (10S). A exatidão é impressionante, particularmente se compararmos a quantidade de trabalho que aqui tivemos com a do Exemplo 6. A integração de Gauss é de considerável importância prática. Sempre que o integrando f for dado por uma fórmula (e não simplesmente por uma tabela de números) ou quando medições experimentais podem ser efetuadas em instantes tj (ou independentemente do que t represente) como mostrado na Tabela 19.7 ou na Ref. [RG1], então a grande acurácia da integração de Gauss se sobrepõe à desvantagem das complicações referentes a tj e Aj (que podem precisar ser armazenados). Além disso, os coeficientes de Gauss Aj são positivos para todo n, diferente- mente de alguns coeficientes de Newton–Cotes para maiores valores de n. É claro que existem freqüentes aplicações com nós igualmente espaçados, casos em que a integração de Gauss não se aplica (ou não representa uma grande vantagem, visto que é preciso primeiro obter-se tj em (11) por interpolação). Como as extremidades –1 e 1 do intervalo de integração em (11) não são zeros de Pn, elas não ocorrem entre os valores t0, • • • , tn, fazendo com que a fórmula de Gauss (11) seja assim chamada de fórmula aberta, para contrastá-la com uma fórmula fechada, onde as extremidades do intervalo de integração são t0 e tn. [Por exemplo, (2) e (7) são fórmulas fechadas.] Tabela 19.7 Integração de Gauss: Nós tj e Coeficientes Aj n Nós tj Coeficientes Aj Grau de precisão 0,57735 02692 1 2 3 0,57735 02692 1 0,77459 66692 0,55555 55556 3 0 0,88888 88889 5 0,77459 66692 0,55555 55556 0,86113 63116 0,34785 48451 4 0,33998 10436 0,65214 51549 7 0,33998 10436 0,65214 51549 0,86113 63116 0,34785 48451 0,90617 98459 0,23692 68851 0,53846 93101 0,47862 86705 5 0 0,56888 88889 9 0,53846 93101 0,47862 86705 0,90617 98459 0,23692 68851 44 Parte E • Análise Numérica 9. Qual é a vantagem das fórmulas de interpolação de Newton em relação às de Lagrange? 10. O que você se lembra sobre os erros na interpolação polinomial? 11. O que é interpolação spline? Qual é a sua vantagem em relação à interpolação polinomial? 12. Liste e compare os métodos numéricos de integração. Quando você os aplicaria? 13. Em que sentido a integração de Gauss é ótima? Explique os deta- lhes. 14. O que significa integração adaptativa? Por que ela é útil? 15. Por que a derivação numérica é geralmente mais delicada que a integração numérica? 16. Escreva –0,35287, 1274,799, –0,00614, 24,9482, 1/3, 85/7 na forma de ponto flutuante com 5S (5 dígitos significativos, apro- priadamente arredondados). 17. Calcule (5,346 – 3,644)/(3,454 – 3,055) conforme dados e depois arredondados passo a passo para 3S, 2S, 1S. (“Passo a passo” significa arredondar os quatro números já arredondados, e não os números dados.) Comente os resultados obtidos. 18. Calcule 0,29731/(4,1232 – 4,0872) com os números conforme dados e depois arredondados passo a passo (ou seja, arredon- dando mais uma vez os números arredondados) para 4S, 3S, 2S. Comente. 19. Resolva x2 – 50x + 1 = 0 por (6) e por (7) na Seção 19.1 usando 5S na computação. Compare e comente. 20. Resolva x2 – 200x + 4 = 0 por (6) e por (7) na Seção 19.1, usando 5S no cálculo. Compare e comente. 21. Considere que os números 4,81 e 12,752 estejam corretamente arredondados ao número de dígitos mostrado. Determine o menor intervalo no qual sua soma (usando os valores verdadeiros ao invés dos arredondados) deve se situar. 22. Responda à questão do Problema 21 para a diferença 4,81 – 12,752. 23. Qual é o erro relativo de na em termos do erro de a? 24. Mostre que o erro relativo de a2 é cerca do dobro do de a. 25. Calcule a solução de x5 = x + 0,2 próximo a x = 0 passando esta equação para a forma x = g(x) e começando de x0 = 0. (Use 6S.) 26. Resolva cos x = x por iteração (6S, x0 = 1), escrevendo-a como x = (0,74x + cos x)/1,74, obtendo x4 = 0,739 085 (valor exato para 6S!). Por que a convergência é tão rápida? 27. Resolva x4 – x3 – 2x – 34 = 0 pelo método de Newton com x0 = 3 e uma precisão de 6S. 28. Resolva cos x – x = 0 pelo método da falsa posição. 29. Calcule f(1,28) a partir de f(1,0) = 3,00000 f(1,2) = 2,98007 f(1,4) = 2,92106 f(1,6) = 2,82534 f(1,8) = 2,69671 f(2,0) = 2,54030 por interpolação linear. Por interpolação quadrática, usando f(1,2), f(1,4), f(1,6). 30. Encontre o spline cúbico para os dados f(–1) = 3 f(1) = 1 f(3) = 23 f(5) = 45 k0 = k3 = 3. 31. Calcule a integral de x3 de 0 a 1 usando a regra do trapézio com n = 5. Que limites de erro são obtidos a partir de (4) da Seção 19.5? Qual é o erro verdadeiro do resultado? Por que esse resultado é maior que o valor exato? 32. Calcule a integral de cos (x2) de 0 a 1 usando a regra de Simpson com 2m = 2 e 2m = 4 e estime o erro usando (10) da Seção 19.5. (Esta é a integral de Fresnel (38) no Apêndice 3.1 com x = 1.) 33. Calcule a integral de cos x de 0 a 12p usando a regra dos três oitavos b a ƒ(x) dx 3 8 h(ƒ0 3ƒ1 3ƒ2 ƒ3) 1 80 (b a)h4ƒ (iv)( t̂ ) e forneça os limites de erro; aqui, a t̂ b e xj = a + (b – a)j/3, j = 0, • • • ,3. R E S U M O D O C A P Í T U L O 1 9 Métodos Numéricos em Geral Neste capítulo, discutimos os conceitos relevantes que perpassam a área de cálculo numérico como um todo, bem como os métodos de natureza geral, em contrapartida aos métodos utilizados em álgebra linear (Capítulo 20) e em equações diferenciais (Capítulo 21). Em cálculos científicos, utilizamos a representação de números em ponto flutuante (Seção 19.1); a representação em ponto fixo é menos adequada na maioria dos casos. Pelos métodos numéricos, obtemos valores aproximados a de quantidades. O erro e de a é (1) e = a – a (Seção 19.1) onde a é o valor exato. O erro relativo de a é e/a. Os erros são causados por arredondamentos, imprecisão de valores medidos, truncamentos (isto é, substituição de integrais por somas, de séries por somas parciais), e assim por diante. Dizemos que um algoritmo é numericamente estável se pequenas alterações nos dados iniciais provocam apenas pequenas alterações correspondentes nos resultados finais. Algoritmos instáveis são em geral inúteis, pois seus erros podem se tornar tão grandes a ponto de deixar os resultados bastante imprecisos. Não se deve confundir a instabilidade Capítulo 19: Métodos Numéricos em Geral 45 numérica dos algoritmos com a instabilidade matemática dos problemas (“problemas mal-condicionados”, Seção 19.2). A iteração de pontos fixos é um método utilizado para resolver equações f(x) = 0, no qual primeiro a equação é transformada algebricamente para x = g(x), faz-se uma suposição inicial x0 para a solução e depois as aproximações x1, x2,• • •, são sucessivamente calculadas por iteração, a partir de (veja a Seção 19.2) (2) xn 1 g(xn) (n 0, 1, • • •). O método de Newton para a solução de equações f(x) = 0 é uma iteração (3) xn 1 xn ƒ(xn) ƒ (xn) (Seção 19.2). Aqui, xn+1 é a interseção entre o eixo x e a tangente à curva y = f(x) no ponto xn. Este método é de segunda ordem (Teorema 2 da Seção 19.2). Se, em (3), substituirmos f por um quociente de diferenças (geometricamente: se substituirmos a tangente por uma secante), obteremos o método da secante; veja (10) na Seção 19.2. Sobre o método da bissecção (que converge lentamente) e o método da falsa posição, veja em Problemas Propostos 19.2. A interpolação polinomial corresponde à determinação de um polinômio pn(x) tal que pn(xj) = fj, onde j = 0,• • •, n e (x0, f0),• • •, (xn, fn) são valores observados ou medidos, ou valores de uma função etc. Dizemos que pn(x) é um polinômio de interpolação. Para dados fornecidos, pn(x) de grau n (ou menos) é único. Entretanto, pode-se escrevê-lo de diferentes formas, notavelmente na forma de Lagrange (4) da Seção 19.3, ou na forma de diferenças divididas de Newton (10) da Seção 19.3, a qual requer um número menor de operações. Para pontos regularmente espaçados x0, x1 = x0 + h,• • •, xn = x0 + nh, esta última passa a ser a fórmula das diferenças progressivas de Newton (fórmula (14) na Seção 19.3) (4) ƒ(x) pn(x) ƒ0 rƒ0 • • • r(r 1) • • • (r n 1) n! nƒ0 onde r = (x – x0)/h, as diferenças progressivas são ∆fj = fj+1 – fj e kƒj k1ƒj 1 k1ƒj (k 2, 3, • • •). Similar a esta é a fórmula de interpolação de diferenças regressivas de Newton (fórmula (18) na Seção 19.3). Os polinômios de interpolação podem se tornar numericamente instáveis à medida que n aumenta, de modo que, ao invés de se interpolar e obter-se uma aproximação por meio de um único polinômio de grau elevado, é preferível utilizar-se um spline cúbico g(x), ou seja, uma função de interpolação dupla e continuamente derivável [assim, g(xj) = fj], que, em cada subintervalo xj x xj+1, consiste em um polinômio cúbico qj(x); veja a Seção 19.4. A regra de Simpson para a integração numérica é [veja (7) na Seção 19.5] (5) b a ƒ(x) dx h 3 (ƒ0 4ƒ1 2ƒ2 4ƒ3 • • • 2ƒ2m2 4ƒ2m1 ƒ2m) com nós igualmente espaçados xj = x0 + jh, j = 1,• • •, 2m, h = (b – a)/(2m), e fj = f(xj). Esta regra é simples, mas precisa o suficiente em diversas aplicações. Seu grau de precisão é GP = 3, pois o erro (8) na Seção 19.5 envolve h4. Uma estimativa de erro mais prática é (10) na Seção 19.5: eh/2 1 15 (Jh/2 Jh), obtida primeiramente calculando-se com um passo h, então com um passo h/2, e então tomando-se 1/15 da diferença entre os resultados. A regra de Simpson é a mais importante das fórmulas de Newton–Cotes, que são obtidas integrando-se os polinômios de interpolação de Lagrange, com os polinômios lineares sendo integrados pela regra do trapézio (2) da Seção 19.5, os polinômios quadráticos pela regra de Simpson, os cúbicos pela regra dos três oitavos (veja em Problemas de Revisão do Capítulo 19) etc. A integração adaptativa (Seção 19.5, Exemplo 6) é uma integração que ajusta (“adapta”) o passo (automaticamente) à variabilidade de f(x). A integração de Romberg (Projeto de Equipe 26, em Problemas Propostos 19.5) parte da regra do trapézio (2) na Seção 19.5, com h, h/2, h/4 etc. e aperfeiçoa os resultados por meio de uma soma sistemática das estimativas de erro. A integração de Gauss (11) na Seção 19.5 é importante devido à sua grande precisão (GP = 2n – 1, comparada com a de Newton–Cotes GP = n – 1 ou n). Tal precisão é alcançada por meio de uma escolha ótima dos nós, que não são igualmente espaçados; veja a Tabela 19.7 na Seção 19.5. A derivação numérica é discutida no final da Seção 19.5. (Sua principal aplicação (em equações diferenciais) é apresentada no Capítulo 21.) C A P Í T U L O 2 0 Métodos Numéricos de Álgebra Linear Neste capítulo, explicaremos alguns dos mais importantes métodos numéricos para a solução de sistemas de equações lineares (Seções 20.1–20.4), para o ajuste com linhas retas ou parábolas (Seção 20.5) e para os problemas com autovalores de matrizes (Seções 20.6–20.9). Esses métodos são de considerável impor- tância prática, pois muitos problemas de engenharia, estatística e outras áreas têm ligação com modelos matemáticos cuja solução requer métodos numéricos de álgebra linear. COMENTÁRIO. Esse capítulo é independente do Capítulo 19 e pode ser estudado imediatamente após os Capítulos 7 ou 8. Pré-requisitos: Seções 7.1, 7.2 e 8.1. Seções que podem ser omitidas em um curso menor: 20.4, 20.5 e 20.9 Referências e Respostas dos Problemas: Parte E do Apêndice 1 e Apêndice 2 20.1 Sistemas Lineares: Eliminação de Gauss Um sistema linear de n equações com n incógnitas x1, • • • , xn consiste em um conjunto de equações E1, • • • , En da forma (1) E1: E2: En: a11x1 • • • a1nxn b1 a21x1 • • • a2nxn b2 • • • • • • • • • • • • • • • • • • • • an1x1 • • • annxn bn onde os coeficientes ajk e bj são números dados. O sistema é chamado de homogêneo se todos os bj são nulos; caso contrário, ele é chamado de não-homogêneo. Usando a multiplicação de matrizes (Seção 7.2), podemos escrever (1) na forma de uma simples equação vetorial (2) Ax b onde a matriz dos coeficientes A = [ajk] é a matriz n n A a11 a21 • an1 a12 a22 • an2 • • • • • • • • • • • • a1n a2n • ann , e x x1 • • • xn e e b b1 • • • bn são vetores-coluna. A seguinte matriz A é chamada de matriz ampliada do sistema (1): A [A b] a11 a21 • an1 • • • • • • • • • • • • a1n a2n • ann b1 b2 • bn . Capítulo 20: Métodos Numéricos de Álgebra Linear 49 a22 5 1 2 2 4, a23 2 1 2 8 2, a24 8 1 2 26 5, e então m31 = 0/6 = 0, de modo que a terceira equação 8x2 + 2x3 = –7 não se alterou no Passo 1. No Passo 2 (k = 2), tivemos 8 como o maior coeficiente na Coluna 2, logo j = 3. Trocamos as posições das equações 2 e 3, calculamos m32 = –4/8 = –12 na linha 4, e a33 = –2 – 1 2 2 = –3, a34 = –5 – 12 ·(–7) = – 3 2 . Isso resultou na forma triangular usada na retrosubstituição. Se, no Passo k, akk = 0, é preciso pivotar. Se akk for pequeno, é recomendável pivotar por causa do aumento dos erros de arredondamento, que pode afetar seriamente a exatidão, ou chegar mesmo a produzir resultados sem sentido. EXEMPLO 3 Dificuldade com Pivôs de Pequeno Valor A solução do sistema 0,0004x1 1,402x2 1,406 0,4003x1 1,502x2 2,501 é x1 = 10, x2 = 1. Resolvemos esse sistema pela eliminação de Gauss, usando uma flutuação pontual aritmética de quatro dígitos. (4D é por razões de simplificação. Faça um exemplo com uma aritmética 8D que mostre o mesmo.) (a) Tomando-se a primeira das equações dadas como a equação-pivô, temos que multiplicar essa equação por m = 0,4003/0,0004 = 1001 e subtrair o resultado da segunda equação, obtendo –1405x2 = –1404. Tabela 20.1 Eliminação de Gauss ALGORITMO DE GAUSS (Ã = [ajk] = [A b]) Este algoritmo calcula uma solução única x = [xj] do sistema (1) ou indica que (1) não possui solução única. ENTRADA: Matriz aumentada n (n + 1) Ã = [ajk], onde aj,n+1 = bj SAÍDA: Solução x = [xj] de (1) ou mensagem dizendo que o sistema (1) não possui solução única Para k = 1, • • • , n – 1, faça: 1 Se ajk = 0 para todo j k então SAÍDA “Não existe uma solução única.” Pare [Procedimento concluído sem sucesso; A é singular] 2 Senão, trocar os conteúdos das linhas j e k de Ã com j sendo o menor j k tal que ajk seja máximo na coluna k. 3 Para j = k + 1, • • • , n, faça: mjk: ajk akk 4 Para p = k + 1, • • • , n + 1, faça: ajp: = ajp – mjkakp Fim Fim Fim 5 Se ann = 0 então SAÍDA “Não existe solução.” Pare Senão 6 xn an,n1 ann [Início da retrosubstituição] Para i = n – 1, • • • , 1, faça: 7 xi 1 aii ai,n1 n ji1 aijxj Fim SAÍDA x = [xj]. Pare Fim GAUSS 50 Parte E • Análise Numérica Logo, x2 = –1404/(–1405) = 0,9993, e da primeira equação, em vez de x1 = 10, temos x1 1 0,0004 (1,406 1,402 0,9993) 0,005 0,0004 12,5. Essa falha ocorre porque a11 é pequeno comparado com a12, de modo que um pequeno erro de arredondamento em x2 conduz a um grande erro em x1. (b) Tomando-se a segunda das equações dadas como a equação-pivô, temos que multiplicar essa equação por 0,0004/0,4003 = 0,000 9993 e subtrair o resultado da primeira equação, obtendo 1,404x2 = 1,404. Logo, x2 = 1 e, da equação-pivô, x1 = 10. Esse cálculo bem-sucedido ocorreu porque a21 não é muito pequeno comparado com a22, de modo que um pequeno erro de arredondamento em x2 não deve conduzir a um grande erro em x1. Com efeito, por exemplo, se tivéssemos o valor x2 = 1,002, ainda chegaríamos pela equação-pivô ao bom valor de x1 = (2,501 + 1,505)/0,4003 = 10,01. A estimativa dos erros para a eliminação de Gauss é discutida na Ref. [E5] no Apêndice 1. O escalonamento por linhas corresponde à multiplicação de cada linha j por um fator adequado de escala sj. Isso é feito em conexão com a pivotação parcial, para se obterem soluções mais acuradas. A despeito da grande quantidade de pesquisas (ver as Refs. [E9] e [E24] no Apêndice 1) e da proposição de diversos princípios, o escalonamento ainda não é uma técnica bem compreendida. Como uma possibilidade, pode-se escalonar somente para a escolha do pivô (e não nos cálculos, a fim de evitar arredondamentos adicionais) e tomar como primeiro pivô a entrada aj1 para a qual aj1/Aj assume o maior valor; aqui Aj é uma entrada de maior valor absoluto na linha j. Algo similar ocorre aos passos adicionais da eliminação de Gauss. Por exemplo, para o sistema 4,0000x1 14020x2 0,4003x1 1,502x2 14060 2,501 poderíamos adotar 4 como pivô, mas dividindo a primeira equação por 104 chegamos no sistema do Exemplo 3, para o qual a segunda equação é uma melhor equação-pivô. Contagem das Operações De modo bastante geral, os fatores importantes para julgar a qualidade de um método numérico são A quantidade de armazenamento A quantidade de tempo ( número de operações) O efeito dos erros de arredondamento. Para a eliminação de Gauss, a contagem de operações para uma matriz cheia (uma matriz com relativamente muitas entradas diferentes de zero) ocorre do seguinte modo. No Passo k, eliminamos xk de n – k equações. Isso exige n – k divisões no cálculo de mjk (linha 3) e (n – k)(n – k + 1) multiplicações, bem como subtrações (ambos na linha 4). Uma vez que executamos n – 1 passos, k vai de 1 a n – 1 e, portanto, o número total de operações nessa etapa da eliminação para frente é ƒ(n) n 1 k 1 (n k) 2 n 1 k 1 (n k)(n k 1) n 1 s 1 s 2 n 1 s 1 s(s 1) 1_2(n 1)n 2_ 3(n 2 1)n 2_3n 3 (escrevemos n – k = s) onde 2n3/3 é obtido descartando-se as potências menores de n. Vemos que f(n) cresce de forma aproximadamente proporcional a n3. Dizemos que f(n) é de ordem n3 e escrevemos ƒ(n) O(n3) onde O significa ordem. A definição geral de O é a seguinte. Escrevemos ƒ(n) O(h(n)) se o quociente f(n)/h(n) permanece limitado (não segue para o infinito) à medida que n → ∞. Em nosso presente caso, h(n) = n3 e, com efeito, f(n)/n3 → 2/3, pois os termos omitidos e divididos por n3 vão para zero à medida que n → ∞. Capítulo 20: Métodos Numéricos de Álgebra Linear 51 Na retrosubstituição de xi fazemos n – i multiplicações e subtrações, além de 1 divisão. Logo, o número de operações na retrosubstituição é b(n) 2 n i 1 (n i) n 2 n s 1 s n n(n 1) n n2 2n O(n2). Vemos que o crescimento é mais devagar que o número de operações na eliminação para a frente do algoritmo de Gauss, de modo que é insignificante para grandes sistemas, pois é menor por um fator n, aproximadamente. Por exemplo, se uma operação consome 10–9 segundos, então os tempos necessários são: Algoritmo n 1.000 n 10.000 Eliminação 0,7 s 11 min Retrosubstituição 0,001 s 0,1 s P R O B L E M A S P R O P O S T O S 2 0 . 1 Para as aplicações de sistemas lineares, veja as Seções 7.1 e 8.2. 1–3 INTERPRETAÇÃO GEOMÉTRICA Resolva graficamente e explique geometricamente. 1. 4x1 3x1 x2 5x2 4,3 33,7 2. 1,820x1 1,183x2 0 12,74x1 8,281x2 0 3. 7,2x1 21,6x1 3,5x2 10,5x2 16,0 48,5 4–14 ELIMINAÇÃO DE GAUSS Resolva os seguintes sistemas lineares pela eliminação de Gauss, com pivotação parcial se necessário (mas sem escalonamento). Mostre os passos intermediários. Cheque o resultado por substituição. Se nenhu- ma solução ou mais de uma solução existir, explique a razão. 4. 6x1 4x1 x2 2x2 3 6 5. 2x1 8x2 6x1 2x2 4 14 6. 7. 8. 9. 20 30 50 10x2 2x3 15x2 3x3 25x2 5x3 4x1 x1 2 3 0 x3 8x3 26x3 3x2 4x2 6x2 5x1 10x1 137,86 85,88 178,54 13x3 8x3 6x2 8x2 6x1 13x1 30,60 8,50 15,48x2 4,30x2 25,38x1 7,05x1 10. 11. 12. 13. 14. 4,65 4,62 4,35 5,97 6,6 x4 8,4 x4 3,0 x4 3,0 x3 5,0x3 7,6x3 4,4 x2 3,6 x2 6,2 x2 0,4 x1 2,0 x1 x1 1,6 32,0 78,0 20,4 7,2 x4 4,8 x4 4,8 x3 9,6 x3 7,2 x3 3,2 x2 1,6 x2 4,8 x2 6,4 x1 3,2 x1 1,20736 2,34066 0,329193 5x3 6x3 3x2 4x23x1 5x1 3,4x1 6,12x2 2,72x3 0 x1 1,80x2 0,80x3 0 2,7x1 4,86x2 2,16x3 0 11 8 2 3x3 x3 2x3 2x2 x2 10x2 x1 10x1 15. EXPERIMENTOS DE ÁLGEBRA COMPUTACIONAL. Eli- minação de Gauss. Escreva um programa para realizar a elimi- nação de Gauss com pivotação. Aplique-o aos Problemas 11–14. Experimente com sistemas onde os determinantes dos coeficientes possuam módulo pequeno. Investigue também o desempenho de seu programa para sistemas maiores de sua escolha, incluindo sis- temas esparsos. 16. PROJETO DE EQUIPE. Sistemas Lineares e Eliminação de Gauss. (a) Existência e unicidade. Encontre a e b tais que ax1 + x2 = b, x1 + x2 = 3 tenham (i) uma única solução, (ii) infinitas soluções e (iii) nenhuma solução. (b) Eliminação de Gauss e inexistência. Aplique a eliminação de Gauss aos seguintes dois sistemas e compare os cálculos 54 Parte E • Análise Numérica 0 11 1 ou 0 1 1 0 não possuem fatoração LU (tente fazer!). Isso indica que, para obter uma fatoração LU, as trocas de posições das linhas de A (e as correspondentes trocas em b) são necessárias. Método de Cholesky Para uma matriz A simétrica, positiva e definida (de modo que, portanto, A = AT, xTAx > 0 para todo x 0), podemos em (2) igualmente escolher U = LT, de modo que ujk = mkj (mas não podemos impor condições aos elementos da diagonal principal). Por exemplo, (5) A LLT . 7 3 5 1 4 0 2 0 0 0 0 5 0 4 3 2 1 7 14 5 83 2 17 5 4 2 14 O método popular de se resolver Ax = b baseado nessa fatoração A = LLT é chamado de método de Cholesky. Em termos dos elementos de L = [ljk], as fórmulas para a fatoração são (6) l11 a11 lj1 j 2, • • • , n lj j aj j j 1 s 1 l js 2 j 2, • • • , n lpj apj j 1 s 1 ljslps p j 1, • • • , n; j 2.1lj j aj1 l11 Se A for simétrica, porém não definida positiva, esse método ainda pode ser aplicado, mas então ele levará a uma matriz L complexa, de modo que o método deixa de ser prático. EXEMPLO 2 Método de Cholesky Resolva pelo método de Cholesky: 4x1 2x1 14x1 2x2 17x2 5x2 14x3 5x3 83x3 14 101 155. Solução. De (6) ou da forma da fatoração, l31 l32 l33 l21 l22 0 l11 0 0 0 0 l33 0 l22 l32 l11 l21 l31 14 5 83 2 17 5 4 2 14 calculamos, na ordem dada, l11 a11 2 l21 1 l31 7 l22 a22 l21 2 17 1 4 l32 (a32 l31l21) ( 5 7 1) 3 l33 a33 l31 2 l32 2 83 72 ( 3)2 5. 1 4 1 l22 14 2 a31 l11 2 2 a21 l11 Isso está de acordo com (5). Agora, temos que resolver Ly = b, isto é, Capítulo 20: Métodos Numéricos de Álgebra Linear 55 . Solução y . 7 27 5 14 101 155 y1 y2 y3 0 0 5 0 4 3 2 1 7 Como segundo passo, temos que resolver Ux = LTx = y, isto é, . Solução x . 3 6 1 7 27 5 x1 x2 x3 7 3 5 1 4 0 2 0 0 TEOREMA 1 Estabilidade da Fatoração de Cholesky A fatoração LLT de Cholesky é numericamente estável (conforme a definição na Seção 19.1). PROVA Temos ajj lj12 lj22 • • • ljj2 ao elevarmos ao quadrado a terceira fórmula em (6) e resolvê-la para ajj. Logo, para todo ljk (note que ljk = 0 para k > j), obtemos (com a desigualdade sendo trivial) ljk2 lj12 lj22 • • • ljj2 ajj. Ou seja, ljk2 é limitado por um elemento de A, o que significa estabilidade contra o arredondamento. Eliminação de Gauss–Jordan. Inversão de Matrizes Outra variante do método da eliminação de Gauss é a eliminação de Gauss–Jordan, introduzida por W. Jordan em 1920, na qual a retrosubstituição é evitada por operações adicionais que reduzem a matriz à forma diagonal, em vez da forma triangular da eliminação de Gauss. Porém, essa redução da forma triangular de Gauss para a forma diagonal requer um maior número de operações do que a retrosubstituição, de modo que esse método é desvantajoso para a solução de sistemas Ax = b. Mas ele pode ser usado para a inversão de matrizes, onde a situação é a seguinte. A inversa de uma matriz quadrada não-singular A pode, em princípio, ser determinada pela resolução de n sistemas (7) Ax bj (j 1, • • • , n) onde bj é a j-ésima coluna da matriz-identidade n n. Entretanto, é preferível chegar a A–1 por meio de operações sobre a matriz-identidade I da mesma forma que no algoritmo de Gauss–Jordan, que reduz A a I. Um típico exemplo que ilustra esse método é dado na Seção 7.8. P R O B L E M A S P R O P O S T O S 2 0 . 2 1–7 MÉTODO DE DOOLITTLE Mostre a fatoração e resolva pelo método de Doolittle. 1. 3x1 15x1 2x2 11x2 15,2 77,3 2. 2x1 9x2 3x1 5x2 41 31 3. 4x1 6x2 34 8x1 7x2 53 4. 2x1 2x1 x1 x2 2x2 2x2 2x3 x3 2x3 0 0 36 5. 6x1 4x2 3x3 4x1 3x2 2x3 3x1 4x2 2x3 2,0 0,5 2,5 6. x1 0,5x1 1,5x1 x2 3,0x2 3,5x2 2,6x3 3,3x3 10,4x3 9,88 16,54 21,02 7. 3x1 18x1 9x1 9x2 48x2 27x2 6x3 39x3 42x3 2,3 13,6 4,5 8. PROJETO DE EQUIPE. O Método de Crout efetua a fatoração A = LU, onde L é triangular inferior e U é triangular superior com elementos da diagonal ujj = 1, j = 1, • • • , n, (a) Fórmulas. Obtenha fórmulas para o método de Crout similares a (4). 56 Parte E • Análise Numérica 20.3 Sistemas Lineares: Solução por Iteração A eliminação de Gauss e suas variantes vistas nas duas últimas seções incluem-se entre os métodos diretos de resolução de sistemas de equações lineares; trata-se de métodos que fornecem soluções após se fazer uma certa quantidade de cálculos que podem ser especificada com antecedência. Em contraste, em um método indireto ou iterativo, começamos a partir de uma aproximação da solução verdadeira e, se obtido sucesso, conseguimos aproximações cada vez melhores por meio de etapas computacionais cíclicas repetidas tantas vezes quantas forem necessárias para se atingir uma exatidão requerida para o resultado, de modo que essa quantidade de aritmética depende da exatidão requerida e varia de caso para caso. Aplicamos métodos iterativos se a convergência for rápida (se as matrizes possuírem grandes valores nos elementos da diagonal principal, conforme veremos), de modo que assim temos uma economia de operações, em comparação com os métodos diretos. Também utilizamos métodos iterativos se um grande sistema for esparso, ou seja, se ele possuir muitos coeficientes nulos, onde poderia haver desperdício de espaço na armazenagem dos zeros, por exemplo, 9995 zeros por equação em um problema de potencial de 104 equações e 104 incógnitas com usualmente apenas 5 termos não-nulos por equação (veja mais sobre isso na Seção 21.4). (b) Exemplos. Resolva os Problemas 1 e 7 pelo método de Crout. (c) Fatore a seguinte matriz pelos métodos de Doolittle, Crout e Cholesky. 1 4 2 4 25 4 2 4 24 (d) Forneça as fórmulas para a fatoração de matrizes tridiagonais pelo método de Crout. (e) Quando é possível obter a fatoração de Crout fazendo-se uma transposição a partir da fatoração de Doolittle? 9–13 MÉTODO DE CHOLESKY Mostre a fatoração e resolva. 9. 9x1 6x1 12x1 6x2 13x2 11x2 12x3 11x3 26x3 87 118 154 10. 0,04x1 0,12x1 0,64x2 0,32x2 0,12x3 0,32x3 0,56x3 1,4 1,6 5,4 11. 4x1 6x1 8x1 6x2 34x2 52x2 8x3 52x3 129x3 0 80 226 12. x1 x1 3x1 2x1 x2 5x2 5x2 2x2 3x3 5x3 19x3 3x3 2x4 2x4 3x4 21x4 30 70 188 2 13. 4x1 2x1 4x1 2x2 4x3 2x2 3x3 3x2 6x3 2x2 3x3 2x4 3x4 9x4 20 36 60 122 14. PROJETO DE ÁLGEBRA COMPUTACIONAL. Método de Cholesky. (a) Escreva um programa para obter a solução de sis- temas lineares pelo método de Cholesky e aplique-o ao Exemplo 2 do texto, aos Problemas 9–11, e a sistemas de sua escolha. (b) Splines. Aplique a parte de fatoração do programa às seguintes matrizes (do modo como elas ocorrem em (9) da Seção 19.4 (com cj = 1), em associação com splines). 2 1 0 1 4 1 0 1 2 , 2 1 0 0 1 4 1 0 0 1 4 1 0 0 1 2 . 15. (Matrizes Definidas) Consideremos que A e B sejam matrizes n n definidas positivas. Serão –A, AT, A + B, A – B definidas positivas? 16–19 INVERSA Encontre a inversa pelo método de Gauss–Jordan, mostrando os deta- lhes do que fizer. 16. No Problema 4. 17. No Problema 5. 18. No Problema 6. 19. No Problema 7. 20. (Arredondamento) Para a seguinte matriz A, encontre det A. O que acontece se arredondarmos os elementos dados para (a) 5S, (b) 4S, (c) 3S, (d) 2S (e) 1S? Qual a implicação prática de seu trabalho? A 1/3 1/9 4/63 1/4 1 3/28 2 1/7 13/49 Capítulo 20: Métodos Numéricos de Álgebra Linear 59 e, multiplicando por (I + L)–1 à esquerda, (7) x(m1) Cx(m) (I L)1b onde C (I L)1U. A iteração de Gauss–Seidel converge para todo x(0) se e somente se todos os autovalores (Seção 8.1) da “matriz de iteração” C = [cjk] tiverem módulos maiores que 1. (A demonstração disso está na Ref. [E5], mencionada no Apêndice 1.) CUIDADO! Se você desejar obter C, primeiramente divida as linhas de A por ajj para fazer com que a diagonal principal seja 1, • • • , 1. Se o raio espectral de C (= valor máximo dos módulos) for pequeno, então a con- vergência é rápida. Condição Suficiente de Convergência. Uma condição suficiente para a convergência é (8) C 1. Aqui, C é alguma norma matricial, como (9) C n j 1 n k 1 c jk 2 (Norma de Frobenius) ou a maior das somas de cjk numa coluna de C (10) C max k n j 1 cjk (Norma da “soma” das colunas) ou a maior das somas de cjk numa linha de C (11) C max j n k 1 cjk (Norma da “soma” das linhas). Essas são as normas de matrizes mais freqüentemente usadas nos métodos numéricos. Na maioria dos casos, a escolha de uma dessas normas é uma questão de conveniência computacional. Entre- tanto, o seguinte exemplo mostra que algumas vezes uma dessas normas é preferível às demais. EXEMPLO 2 Teste de Convergência da Iteração de Gauss–Seidel Teste se a iteração de Gauss–Seidel converge para o sistema 2x x x y 2y y z 4 z 4 2z 4 escrevendo como x 2 12y 1 2z y 2 12x 1 2z z 2 12x 1 2y. Solução. A decomposição (multiplique a matriz por 1/2 – por quê?) é I L U I . 1/2 1/2 0 1/2 0 0 0 0 0 0 0 0 0 0 1/2 0 1/2 1/2 1/2 1/2 1 1/2 1 1/2 1 1/2 1/2 Isso mostra que C (I L) 1U . 1/2 1/4 3/8 1/2 1/4 1/8 0 0 0 1/2 1/2 0 1/2 0 0 0 0 0 0 0 1 0 1 1/2 1 1/2 1/4 Calculamos a norma de Frobenius de C C 1/2 1/2 0,884 1 50 64 9 64 1 64 1 16 1 16 1 4 1 4 e concluímos de (8) que essa iteração de Gauss–Seidel converge. É interessante que as outras duas normas não permitiriam essa conclusão, como você deveria verificar. Obviamente, isso aponta para o fato de que (8) é suficiente para a convergência, e não apenas necessária. Resíduo. Dado um sistema Ax = b, o resíduo r de x com relação a esse sistema é definido por 60 Parte E • Análise Numérica (12) r b Ax. Logicamente, r = 0 se e somente se x é uma solução. Logo r 0 para uma solução aproximada. Na iteração de Gauss–Seidel, a cada estágio modificamos ou relaxamos uma componente de uma solução aproximada, no intuito de reduzir uma componente de r a zero. Logo, a iteração de Gauss–Seidel pertence a uma classe de métodos freqüentemente chamada de métodos de relaxação. Maiores informações sobre os resíduos são fornecidas na próxima seção. Iteração de Jacobi A iteração de Gauss–Seidel é um método de correções sucessivas, pois, para cada componente, substituímos sucessivamente uma aproximação de uma componente por uma nova aproximação correspondente, tão logo esta última tenha sido calculada. Dizemos que um método iterativo é de correções simultâneas se nenhuma compo- nente da aproximação x(m) é usada até que todas as componentes de x(m) tenham sido calculadas. Um método desse tipo é a iteração de Jacobi, que é similar à iteração de Gauss–Seidel, mas envolve a não-utilização de valores melhorados até que um passo tenha sido concluído e então substituímos x(m) por x(m+1) uma vez, imediatamente antes do próximo passo. Logo, se escrevemos Ax = b (com ajj = 1 como antes!) na forma x = b + (I – A)x, a iteração de Jacobi em notação matricial é (13) x(m1) b (I A)x(m) (ajj 1). Esse método converge para toda escolha de x(0) se e somente se o raio espectral de I – A for menor que 1. Isso vem recentemente ganhando grande interesse prático, devido ao fato de que, em processos paralelos, todas as n equações podem ser resolvidas simultaneamente a cada passo da iteração. Sobre Jacobi, veja a Seção 10.3. Sobre exercícios, veja Problemas Propostos. P R O B L E M A S P R O P O S T O S 2 0 . 3 1. Verifique a declaração feita no final do Exemplo 2. 2. Mostre que, para o sistema do Exemplo 2, a iteração de Jacobi diverge. Sugestão: use autovalores. 3–8 ITERAÇÃO DE GAUSS–SEIDEL Execute 5 passos, começando de x0 = [1 1 1]T e usando 6S na com- putação. Sugestão: certifique-se de ter resolvido cada equação para a variável com o maior coeficiente (por quê?). Mostre os detalhes. 3. 4. 5. 19 2 39 2x3 2x3 8x3 x2 4x2 3x2 5x1 x1 2x1 25,5 0 10,5 7x3 x3 x2 x2 6x2 5x1 x1 61,3 49,1 185,8 6x3 2x3 x3 x2 9x2 2x2 x1 x1 8x1 6. 21 45 33 x3 4x3 x2 4x2 x2 4x1 x1 7. 8. 12,5 18,5 11,5 5x3 2x3 x3 6x2 2x2 4x1 x1 8x1 x3 6 x3 6 10x3 6 x2 10x2 x2 10x1 x1 x1 9. Aplique a iteração de Gauss–Seidel (3 passos) ao sistema do Pro- blema 7, começando por (a) 0, 0, 0, (b) 10, 10, 10. Compare e comente. 10. No Problema 7, calcule C (a) se você resolver a primeira equação para x1, a segunda para x2, a terceira para x3, provando a convergên- cia. (b) se você, equivocadamente, resolver a terceira equação para x1, a primeira para x2, a segunda para x3, provando a divergência. 11. PROJETO DE ÁLGEBRA COMPUTACIONAL. Iteração de Gauss–Seidel. (a) Escreva um programa para a fazer a iteração de Gauss–Seidel. (b) Aplique o programa a A(t)x = b, começando de [0 0 0]T, onde A(t) , b . 2 2 2 t t 1 t 1 t 1 t t Capítulo 20: Métodos Numéricos de Álgebra Linear 61 20.4 Sistemas Lineares: Mau Condicionamento, Normas Não se precisa de muita experiência para observar que alguns sistemas Ax = b são bons, fornecendo soluções exatas mesmo com arredondamentos ou coeficientes inexatos, ao passo que outros são ruins, onde as inexatidões afetam fortemente a solução. Queremos ver o que acontece, e se é ou não possível “confiar” num sistema linear. Vamos primeiramente formular dois conceitos relevantes (os de mau e de bom condicionamento) para desenvolvi- mentos numéricos em geral e então retornar aos sistemas lineares e matrizes. Dizemos que um problema computacional é mal condicionado (ou mal posicionado) se “pequenas” altera- ções nos dados (a entrada) provocam “grandes” alterações na solução (a saída). Por outro lado, dizemos que um problema é bem condicionado (ou bem posicionado) se “pequenas” alterações nos dados provocam apenas “pequenas” alterações na solução. Esses conceitos são qualitativos. Certamente, consideraríamos “grande” um aumento das imprecisões por um fator de 100, porém poderia haver controvérsias quanto à determinação exata da linha que separa o que é “gran- de” do que é “pequeno”, dependendo do tipo de problema e do nosso ponto de vista. Uma precisão dupla pode algumas vezes ajudar, mas se os dados são medidos de forma imprecisa, deveríamos tentar modificar o ajuste matemático do problema, de modo a termos um problema bem condicionado. Retornemos agora aos sistemas lineares. A Fig. 442 explica que o mau condicionamento ocorre se e somente se as duas equações fornecem duas retas aproximadamente paralelas, de modo que seu ponto de interseção (a Para t = 0,2, 0,5, 0,8, 0,9, determine o número de passos necessários para a obtenção da solução exata, usando 6S e o correspondente raio espectral de C. Faça um gráfico do número de passos e do raio espectral como funções de t e comente. (c) Sobre-relaxação sucessiva (SOR*). Mostre que, adicionando e subtraindo x(m) à direita, a fórmula (6) pode ser escrita como x(m1) x(m) b Lx(m1) (U I)x(m) (ajj 1). A possível ocorrência de correções adicionais motiva a introdução de um fator de sobre-relaxação ω > 1 para assim se obter a fór- mula de SOR para Gauss–Seidel (14) x(m1) x(m) v(b Lx(m1) (U I)x(m)) (ajj 1) com a finalidade de fornecer uma convergência mais rápida. Um valor recomendado é v 2/(1 1 r) onde r é o raio espectral de C em (7). Aplique a sobre-relaxação sucessiva à matriz em (b) para t = 0,5 e 0,8 e observe a melhoria da convergência. (Ganhos espetaculares são obtidos com sistemas maiores.) 12–15 ITERAÇÃO DE JACOBI Execute 5 passos, começando de x0 = [1 1 1]T. Compare com a ite- ração de Gauss–Seidel. Com qual dos dois métodos a convergência parece ser mais rápida? (Mostre os detalhes do que fizer.) 12. O sistema do Problema 6 13. O sistema do Problema 5 14. O sistema do Problema 8 15. Mostre a convergência no Problema 14 verificando que I – A, onde A é a matriz do Problema 14 com as linhas divididas pelos cor- respondentes elementos da diagonal principal, tem os autovalores –0,519589 e 0,259795 ± 0,246603i. 16–20 NORMAS Calcule as normas (9), (10), (11) para as seguintes matrizes (quadra- das). Comente as razões das maiores ou menores diferenças entre os três números. 16. A matriz do Problema 3 17. A matriz do Problema 7 18. A matriz do Problema 8 19. 20. 5 0 2 3 1 12 4 7 17 k k 2k k 2k k 2k k k *Abreviatura da expressão original: successive overrelaxation. (N.T.) y x γ y x (a) (b) Fig. 442. Sistema linear de duas equações e duas incógnitas (a) bem condicionado e (b) mal condicionado 64 Parte E • Análise Numérica Tomando nosso melhor c = A possível (ou seja, nosso menor c), temos de (8) (11) Ax A x. Essa é a fórmula de que precisaremos. Para as matrizes n n (veja a Ref. [E17]), a fórmula (9) também implica que (12) AB A B, portanto, An An. Para outras fórmulas usuais de normas, veja as Refs. [E9] e [E17]. Antes de prosseguirmos, façamos um cálculo simples a título de exemplo. EXEMPLO 4 Normas Matriciais Calcule as normas matriciais da matriz A dos coeficientes no Exemplo 1 e de sua inversa A–1, supondo que usamos (a) a norma vetorial l1 e (b) a norma vetorial l∞. Solução. Usamos (4*) da Seção 7.8 para a inversa e então (10) e (11) da Seção 20.3. Portanto, A 0,9999 1,00011,0000 1,0000 , A1 5000,0 5000,5 5000,0 4999,5 . (a) A norma vetorial l1 fornece a norma da “soma” das colunas (10) da Seção 20.3; da coluna (2), obtemos, portanto, A = –1,0001 + –1,0000 = 2,0001. De forma similar, A–1 = 10.000. (b) A norma vetorial l∞ fornece a norma da “soma” das linhas (11) da Seção 20.3; portanto, da linha 1, A = 2, A–1 = 10.000,5. Notamos que A–1 é surpreendentemente grande, fazendo com que o produto A A–1 seja grande (20.001). Veremos a seguir que isso usualmente ocorre nos sistemas mal condicionados. Número Condicional de uma Matriz Estamos agora prontos para introduzir o conceito-chave da nossa discussão sobre mau condicionamento, a saber, o número condicional k(A) de uma matriz quadrada (não-singular) A, definido por (13) k(A) A A1. O papel do número condicional é apresentado pelo seguinte teorema. TEOREMA 1 Número Condicional Um sistema linear de equações Ax = b e sua matriz A cujo número condicional (13) seja pequeno são bem condicionados. Um grande número condicional indica mau condicionamento. PROVA b = Ax e (11) fornecem b A x. Consideremos que b 0 e x 0. Então, dividindo por b x, temos (14) 1 x A b . Multiplicando (2) r = A(x – x) por A–1 no lado esquerdo e invertendo os lados, temos x – x = A–1r. Agora, (11), com A–1 e r ao invés de A e x, fornece x x A1r A1 r. Dividindo por x [note que x 0 por (3b)] e usando (14), finalmente obtemos (15) x x x 1 x A1 r A b A1 r k(A) r b . Logo, se k(A) é pequeno, um pequeno r / b implica um pequeno erro relativo x – x / x, de modo que o sistema é bem condicionado. Entretanto, isso não ocorrerá se k(A) for grande; nesse caso, um pequeno r / b não necessariamente implicará um pequeno erro relativo x – x / x. Capítulo 20: Métodos Numéricos de Álgebra Linear 65 EXEMPLO 5 Números Condicionais. Iteração de Gauss–Seidel A 5 1 1 1 4 2 1 2 4 tem a inversa A1 156 12 2 2 2 19 9 2 9 19 . Como A é simétrica, (10) e (11) na Seção 20.3 fornecem o mesmo número condicional k(A) A A1 7 156 30 3,75. Vemos que um sistema linear Ax = b com esse A é bem condicionado. Por exemplo, se b = [14 0 28]T, o algoritmo de Gauss fornece a solução x = [2 –5 9]T (confirme isso). Como os elementos da diagonal principal de A são relativamente grandes, podemos esperar uma convergência razoavelmente boa para a iteração de Gauss–Seidel. Com efeito, começando de, digamos, x0 = [1 1 1]T, obtemos os 8 primeiros passos (com valores 3D) x1 x2 x3 1,000 1,000 1,000 2,400 1,100 6,950 1,630 3,882 8,534 1,870 4,734 8,900 1,967 4,942 8,979 1,993 4,988 8,996 1,998 4,997 8,999 2,000 5,000 9,000 2,000 5,000 9,000 EXEMPLO 6 Sistema Linear Mal Condicionado Por (10) ou (11) da Seção 20.3, o Exemplo 4 fornece para a matriz do Exemplo 1 o grande número condicional k(A) = 2,0001 10 000 = 2 10 000,5 = 20 0001. Isso confirma o fato de que o sistema é bastante mal condicionado. Algo similar ocorre no Exemplo 2, onde, por (4*) da Seção 7.8 e com um cálculo em 6D, A1 1 0,0002 1,00011,0000 1,0000 1,0001 5000,5 5000,0 5000,0 5000,5 de modo que (10) da Seção 20.3 conduz a um k(A) muito grande, explicando o surpreendente resultado do Exemplo 2, k(A) (1,0001 1,0000)(5000,5 5000,0) 20 002. Na prática, não conheceremos A–1, de modo que, ao calcularmos o número condicional k(A), precisamos estimar A–1. Um método de fazer isso (proposto em 1979) é explicado na Ref. [E9] listada no Apêndice 1. Elementos Inexatos da Matriz. k(A) pode ser usado para estimar o efeito δx de uma imprecisão δA de A (erros de medições dos ajk, por exemplo). Em vez de Ax = b, temos, então, (A dA)(x dx) b. Multiplicando e subtraindo Ax = b em ambos os lados, obtemos Adx dA(x dx) 0. Multiplicando por A–1 à esquerda e passando o segundo termo para a direita, temos dx A1dA(x dx). Aplicando (11) com A–1 e o vetor δA(x + δx) em vez de A e x, obtemos dx A1dA(x dx) A1 dA(x dx). Aplicando (11) à direita, com δA e x – δx em vez de A e x, obtemos dx A1 dA x dx. Agora, pela definição de k(A), A–1 = k(A)/A, de modo que a divisão por x + δx mostra que o erro relativo de x relaciona-se ao de A por meio do número condicional, através da desigualdade 66 Parte E • Análise Numérica (16) dx x dx x dx A1 dA k(A) dA A . Conclusão. Se o sistema for bem condicionado, pequenas inexatidões δA / A podem exercer somente um pequeno efeito sobre a solução. Por outro lado, em caso de mau condicionamento, se δA / A for pequeno, δx / x poderá ser grande. Inexatidão do Lado Direito. Você pode mostrar que, de forma similar, quando A é exata, uma inexatidão δb de b provoca uma inexatidão δx que satisfaz a (17) dx x k(A) db b . Logo, δx / x precisa permanecer relativamente pequeno sempre que k(A) for pequeno. EXEMPLO 7 Inexatidões. Limites (16) e (17) No Exemplo 5, se cada um dos nove elementos de A for medido com uma imprecisão de 0,1, então δA = 9 0,1 e (16) fornece dx x 7,5 3 0,1 7 0,321 portanto, dx 0,321 x 0,321 16 5,14. Você encontrará experimentalmente que a inexatidão real δx é apenas cerca de 30% do limite 5,14. Isso é usual. De forma similar, também no Exemplo 5, se δb = [0,1 0,1 0,1]T, então δb = 0,3 e b = 42, de modo que (17) fornece dx x 7,5 0,3 42 0,0536, logo, dx 0,0536 16 0,857 porém esse limite é novamente muito maior que a inexatidão atual, que vale cerca de 0,15. Comentários Adicionais Sobre os Números Condicionais. As seguintes explicações adicionais podem ser úteis. 1. Não há uma linha precisa separando o “bom condicionamento” do “mal condicionamento”, mas geralmente a situação vai piorando à medida que vamos passando de sistemas com pequenos valores de k(A) para sistemas com grandes valores de k(A). Agora, sempre ocorre que k(A) 1, de modo que valores de 10, 20, ou algo parecido não representam qualquer motivo de preocupação; por outro lado, a ocorrência de, digamos, k(A) = 100 exige cautela, e sistemas como os dos Exemplos 1 e 2 são extremamente mal condicionados. 2. Se k(A) for grande (ou pequeno) em uma norma, ele será grande (ou pequeno, respectivamente) em qualquer outra norma. Veja o Exemplo 5. 3. A literatura sobre o mau condicionamento é extensa. Para uma introdução a esse assunto, veja [E9]. Isto encerra nossa discussão sobre os métodos numéricos para a solução de sistemas lineares. Na próxima seção, consideraremos o ajuste de curvas, uma importante área que permite a obtenção de soluções de sistemas lineares. P R O B L E M A S P R O P O S T O S 2 0 . 4 1–8 NORMAS VETORIAIS Compute (5), (6), (7). Calcule um vetor unitário correspondente (vetor de norma 1) com relação à norma l∞. 1. [1 –6 5] 2. [0,4 –1,2 0 8,0] 3. [–4 4 3 –3] 4. [0 0 1 0 0] 5. [0,3 –0,1 0,5 1,0] 6. [16 21 54 –119] 7. [1 1 1 1 1 1] 8. [3 0 0 –3 0] 9. Mostre que x∞ x2 x1. 10–15 NORMAS MATRICIAIS, NÚMEROS CONDICIONAIS Calcule a norma matricial e o número condicional correspondente à norma vetorial l1. 10. 11. 12. 13. 14. 5,25 4,2 3,5 3 7 5,25 4,2 3,5 10,5 7 5,25 4,2 21 10,5 7 5,25 100 0 0 0 100 0 0 0 0,01 3 3 3 0 7 10 5 7 4 2 3 1 Capítulo 20: Métodos Numéricos de Álgebra Linear 69 EXEMPLO 1 Linha Reta Usando o método dos mínimos quadrados, ajuste uma linha reta aos quatro pontos dados na fórmula (1). Solução. Obtemos n 4, xj 0,1, x j 2 3,43, yj 3,907, xjyj 2,3839. Logo, as equações normais são 4a 0,10b 3,9070. 0,1a 3,43b 2,3839. A solução (arredondada para 4D) é a = 0,9601, b = 0,6670, e obtemos a linha reta (Fig. 443) y 0,9601 0,6670x. Ajuste de Curvas por Polinômios de Grau m Nosso método de ajuste de curvas pode ser generalizado de um polinômio y = a + bx para um polinômio de grau m (5) p(x) b0 b1x • • • bmxm onde m n – 1. Então, q assume a forma q n j1 (yj p(xj))2 e depende de m + 1 parâmetros b0, • • • , bm. Em vez de (3), temos então m + 1 condições (6) q b0 0, • • • , q bm 0 que fornecem um sistema de m + 1 equações normais. No caso de um polinômio quadrático (7) p(x) b0 b1x b2x2 as equações normais são (somando-se de 1 a n) (8) b0n b0 xj b0 x j2 b1 xj b2 x j2 yj b1 x j2 b2 x j3 xjyj b1 x j3 b2 x j4 x j2 yj. Deixamos a obtenção de (8) a cargo do leitor. EXEMPLO 2 Parábola Quadrática por Mínimos Quadrados Ajuste uma parábola que passe através dos dados (0, 5), (2, 4), (4, 1), (6, 6), (8, 7). Solução. Para as equações normais, precisamos de n 5, xj 20, x j2 120, x j3 800, x j4 5664 yj 23, xjyj 104, x j2yj 696. Logo, essas equações são 5b0 20b0 120b0 20b1 120b1 800b1 120b2 800b2 5664b2 23 104 696. Resolvendo-as, obtemos a parábola quadrática de mínimos quadrados (Fig. 445) y 5,11429 1,41429x 0,21429x2. Para um polinômio geral (5), as equações normais formam um sistema linear de equações nas incógnitas b0, • • • , bm. Quando sua matriz M é não-singular, podemos resolver o sistema pelo método de Cholesky (Seção 20.2), pois M será então definida positiva (e simétrica). Quando as equações são de modo aproximado linearmente depen- dentes, as equações normais podem tornar-se mal condicionadas, devendo ser substituídas por outros métodos; veja [E5] da Seção 5.7, listada no Apêndice 1. O método dos mínimos quadrados também desempenha um papel em estatística (veja a Seção 25.9). 70 Parte E • Análise Numérica 1–6 AJUSTANDO UMA LINHA RETA Usando o método dos mínimos quadrados, ajuste uma linha reta para os pontos dados. Mostre os detalhes. Verifique seu resultado esboçan- do os pontos e a linha. Julgue a eficácia do ajuste. 1. (2, 0), (3, 4), (4, 10), (5, 16) 2. Como a linha do Problema 1 se altera se adicionamos um ponto muito acima dela, digamos (3, 20)? 3. (2,5, 8,0), (5,0, 6,9), (7,5, 6,2), (10,0, 5,0) 4. (Lei de Ohm U = Ri) Pela linha dos mínimos quadrados, estime a resistência R que ajusta (i, U) = (2,0, 104), (4,0, 206), (6,0, 314), (10,0, 530). 5. (Velocidade média) Estime a velocidade média vav de um carro viajando de acordo com s = v t [km] (s = distância percorrida, t [h] = tempo) de (t, s) = (9, 140), (10, 220), (11, 310), (12, 410). 6. (Lei de Hooke F = ks) Estime a constante elástica k de uma mola a partir da força F [lb] e da elongação s [cm], onde (F, s) = (1, 0,50), (2, 1,02), (4, 1,99), (6, 3,01), (10, 4,98), (20, 10,03). 7. Demonstre as equações normais (8). 8–10 AJUSTANDO UMA PARÁBOLA QUADRÁTICA Usando mínimos quadrados, ajuste uma parábola (7) para os pontos dados (x, y). Verifique o resultado por meio de um esboço. 8. (–1, 3), (0, 0), (1, 2), (2, 8) 9. (0, 4), (2, 2), (4, –1), (6, –5) 10. Tempo x [h] do funcionário encarregado 1 2 3 4 5 Tempo [s] de reação do funcionário 1,50 1,28 1,40 1,85 2,20 11. Ajuste (2) e (7) usando mínimos quadrados para (–1,0, 5,4), (–0,5, 4,1), (0, 3,9), (0,5, 4,8), (1,0, 6,3), (1,5, 9,3). Num mesmo gráfico, represente os dados e as curvas e comente. 12. (Parábola cúbica) Usando o método dos mínimos quadrados, obtenha a fórmula para as equações normais de uma parábola cúbica. 13. Usando o método dos mínimos quadrados, ajuste as curvas (2) e (7) e uma parábola cúbica para (–2, –35), (–1, –9), (0, –1), (1, –1), (2, 17), (3, 63). Represente num mesmo gráfico as três curvas e os pontos. Comente a eficácia do ajuste. 14. PROJETO DE ÁLGEBRA COMPUTACIONAL. Mínimos Quadrados. Escreva programas para calcular e resolver as equa- ções normais (4) e (8). Aplique os programas aos Problemas 3, 5, 9, 11. Se o seu aplicativo tiver um comando específico para ajuste (como têm o Maple e o Mathematica), compare seus resultados com os obtidos pelo seu programa. y x0 2 4 6 8 2 4 6 8 Fig. 445. Parábola dos mínimos quadrados do Exemplo 2 P R O B L E M A S P R O P O S T O S 2 0 . 5 15. EXPERIMENTO DE ÁLGEBRA COMPUTACIONAL. Míni- mos Quadrados versus Interpolação. Para os dados fornecidos e também para dados que você escolher, encontre o polinômio de interpolação e as aproximações por mínimos quadrados (linear, quadrática etc.). Compare e comente. (a) (–2, 0), (–1, 0), (0, 1), (1, 0), (2, 0) (b) (–4, 0), (–3, 0), (–2, 0), (–1, 0), (0, 1), (1, 0), (2, 0), (3, 0), (4, 0) (c) Escolha cinco pontos sobre uma linha reta, por exemplo, (0, 0), (1, 1), • • • , (4, 4). Mova um dos pontos uma unidade para cima e encontre o polinômio quadrático pelo método dos mínimos quadrados. Faça isso para cada ponto. Represente num mesmo gráfico os cinco polinômios. Qual dos cinco movimentos tem o maior efeito? 16. PROJETO DE EQUIPE. A aproximação por mínimos qua- drados de uma função f(x) num intervalo a x b por uma função Fm(x) a0y0(x) a1y1(x) • • • amym(x) onde y0(x), • • • , ym(x) são funções dadas, requer a determinação dos coeficientes a0, • • • , am tais que (9) b a [ƒ(x) Fm(x)]2 dx torna-se mínimo. Essa integral é denotada por f – Fm2, e f – Fm é chamada de norma L2 de f – Fm (com L em homenagem a Lebes- gue2). Uma condição necessária para esse mínimo é dada por ∂f – Fm2/∂aj = 0, j = 0, • • • , m [o análogo de (6)]. (a) Mostre que isso conduz a m + 1 equações normais (j = 0, • • • , m) m k0 hjkak bj onde (10) hjk b a yj(x)yk(x) dx, bj b a ƒ(x)yj(x) dx. (b) Polinômio. Que forma assume (10) se Fm(x) a0 a1x • • • amxm ? Qual é a matriz dos coeficientes de (10) nesse caso quando o intervalo é 0 x 1? (c) Funções ortogonais. Quais são as soluções de (10) se y0(x), • • • , ym(x) são ortogonais no intervalo a x b? (Para a definição, veja a Seção 5.7. Veja também a Seção 5.8.) 2HENRI LEBESGUE (1875–1941), grande matemático francês, que criou a moderna teoria de medidas e integração em sua famosa tese de doutoramento de 1902. Capítulo 20: Métodos Numéricos de Álgebra Linear 71 20.6 Problemas de Matrizes de Autovalores: Introdução Nas seções restantes deste capítulo, discutiremos algumas das mais importantes idéias e métodos numéricos para problemas de matrizes de autovalores. Essa extensa parte da álgebra linear numérica tem grande importância prática e sobre ela há muita pesquisa sendo desenvolvida, com centenas, senão milhares, de artigos publicados em vários periódicos especializados (veja as Refs. em [E8], [E9], [E11], [E29]). Começaremos com os conceitos e os resultados gerais de que precisaremos para explicar e aplicar os métodos numéricos aos problemas de autovalores. (Para modelos típicos de problemas de autovalores, veja o Capítulo 8.) Um autovalor ou valor característico (ou ainda raiz latente) de uma dada matriz n n A = [ajk] é um número real ou complexo λ tal que a equação vetorial (1) Ax lx possui uma solução não-trivial, ou seja, uma solução x 0, que é então chamada de autovetor ou vetor carac- terístico da A correspondente a esse autovalor λ. O conjunto de todos os autovalores de A é chamado de espectro de A. A Equação 1 pode ser escrita como (2) (A lI)x 0 onde I é a matriz-identidade n n. Esse sistema homogêneo possui uma solução não-trivial se e somente se o determinante característico det(A – λI) é 0 (veja o Teorema 2 da Seção 7.5). Isso nos leva ao (veja a Seção 8.1). TEOREMA 1 Autovalores Os autovalores de A são as soluções λ da equação característica (3) det (A lI) a11 l a21 • an1 a12 a22 l • an2 • • • • • • • • • • • • a1n a2n • ann l 0. Desenvolvendo o determinante característico, obtemos o polinômio característico de A, que é de grau n em λ. Logo, A tem pelo menos um e no máximo n diferentes números de autovalores. Se A é real, o mesmo ocorre com os coeficientes do polinômio característico. Da álgebra comum, decorre que, então, as raízes (os autovalores de A) são reais ou conjugados complexos aos pares. Usualmente representaremos os autovalores de A por l1, l2, • • • , ln porém admitindo que alguns deles (ou todos) podem ser iguais. A soma desses n autovalores é igual à soma dos elementos da diagonal principal de A, chamada de traço de A; portanto, (4) traço A n j1 aj j n k1 lk. Além disso, o produto dos autovalores é igual ao determinante de A, (5) det A l1l2 • • • ln. Ambas as fórmulas decorrem da fatoração do polinômio característico, que denotamos por f(λ), ƒ(l) (1)n(l l1)(l l2) • • • (l ln). 74 Parte E • Análise Numérica D1: Centro 0, raio 1, D2: Centro 5, raio 1,5, D3: Centro 1, raio 1,5. Os centros são os elementos da diagonal principal de A. Eles seriam os autovalores de A se A fosse diagonal. Podemos tomar esses valores como aproximações grosseiras dos autovalores desconhecidos (valores 3D) λ1 = –0,209, λ2 = 5,305, λ3 = 0,904 (verifique isso); então, os raios dos discos correspondem aos limites dos erros. Como A é simétrica, do Teorema 5 da Seção 20.6 segue que o espectro de A precisa de fato se situar nos intervalos [–1, 2,5] e [3,5 , 6,5]. É interessante que, aqui, os discos de Gerschgorin formam dois conjuntos disjuntos, a saber, D1 D3, que contêm dois autovalores, e D2, que contém um autovalor. Isso é usual, conforme mostra o seguinte teorema. 0 D 1 D 2 D 3 1 5 y x Fig. 446. Discos de Gerschgorin do Exemplo 1 TEOREMA 2 Extensão do Teorema de Gerschgorin Se, numa dada matriz A, p discos de Gerschgorin formam um conjunto S que é disjunto de n – p outros discos, então S contém precisamente p autovalores de A (cada um com sua multiplicidade algébrica, con- forme definida na Seção 20.6). Idéia da Prova. Consideremos A = B + C, onde B é a matriz diagonal com elementos ajj, e apliquemos o Teo- rema 1 a At = B + tC com t real crescendo de 0 a 1. EXEMPLO 2 Outras Aplicações do Teorema de Gerschgorin. Similaridade Suponha que diagonalizamos a matriz por algum método numérico que nos deixou com alguns elementos de fora da diagonal com um tamanho de 10–5, como, por exemplo, A 2 105 105 105 2 105 105 105 4 . O que podemos concluir sobre os desvios dos autovalores em relação aos elementos da diagonal principal? Solução. Pelo Teorema 2, um autovalor deve se situar no disco de raio 2 10–5 centrado em 4 e dois autovalores (ou um autovalor de multiplicidade algébrica 2) no disco de raio 2 10–5 centrado em 2. Com efeito, como a matriz é simétrica, esses autovalores necessariamente se situam nas interseções desses discos com o eixo real, segundo o Teorema 5 da Seção 20.6. Mostremos de que modo um disco isolado pode sempre ser reduzido em tamanho por uma transformação de similaridade. A matriz B T1AT 1 0 0 0 1 0 0 0 105 2 105 105 105 2 105 105 105 4 1 0 0 0 1 0 0 0 105 2 105 1010 105 2 1010 1 1 4 é similar a A. Logo, pelo Teorema 2 da Seção 20.6, ela possui os mesmos autovalores de A. Da linha 3 temos o menor disco de raio 2 10–10. Note que os outros discos ficaram maiores, aproximadamente por um fator 105. E, ao escolhermos T, temos que observar que os novos discos não se sobrepõem ao disco cujo tamanho desejamos diminuir. Para mais fatos interessantes, veja o livro [E28], recentemente publicado. Por definição, uma matriz diagonalmente dominante A = [ajk] é uma matriz n n tal que (3) ajj kj ajk j 1, • • • , n Capítulo 20: Métodos Numéricos de Álgebra Linear 75 onde o somatório engloba todos os elementos de fora da diagonal na linha j. Dizemos que a matriz é estrita e diago- nalmente dominante se > em (3) para todo j. Use o Teorema 1 para provar a seguinte propriedade básica. TEOREMA 3 Dominância Estrita e Diagonal Matrizes estrita e diagonalmente dominantes são não-singulares. Outros Teoremas de Inclusão Um teorema de inclusão é aquele que especifica um conjunto que contém ao menos um autovalor de uma dada matriz. Portanto, os Teoremas 1 e 2 são de inclusão e chegam mesmo a incluir a totalidade do espectro. Discutiremos agora alguns famosos teoremas que produzem outras inclusões de autovalores. Enunciaremos sem demonstração os dois primeiros deles (pois prová-los excederia o nível deste livro). TEOREMA 4 Teorema de Schur3 Consideremos A = [ajk] uma matriz n n. Então, para cada um de seus autovalores λ1, • • • , λn, (4) lm2 n i1 li2 n j1 n k1 ajk2 (Desigualdade de Schur). Em (4), o segundo sinal de igualdade se verifica se e somente se A seja tal que (5) A TA AA T. As matrizes que satisfazem (5) são chamadas de matrizes normais. Não é difícil ver que as matrizes hermitia- nas, anti-hermitianas e unitárias são normais, o mesmo ocorrendo com as matrizes simétricas, anti-simétricas e ortogonais reais. EXEMPLO 3 Limites para os Autovalores Obtidos pela Desigualdade de Schur Para a matriz A 26 2 4 2 21 2 2 4 28 Obtemos, pela desigualdade de Schur, l 1949 44,1475. Você pode verificar que os autovalores são 30, 25, e 20. Então, 302 + 252 + 202 = 1925 < 1949; de fato, A não é normal. Os teoremas precedentes são válidos para toda matriz quadrada real ou complexa. Outros teoremas aplicam-se somente a classes especiais de matrizes. O seguinte é famoso. TEOREMA 5 Teorema de Perron4 Consideremos que A seja uma matriz real n n cujos elementos são todos positivos. Então, A tem um autovalor real positivo λ = r de multiplicidade 1. O autovetor correspondente pode ser escolhido com todas as componentes positivas. (Os outros autovalores são menores que ρ em valor absoluto.) Para a demonstração, veja a Ref. [B3], vol. II. O teorema também se aplica às matrizes com termos reais não- negativos (“Teorema de Perron–Frobenius”4), desde que A seja irredutível, isto é, não possa ser colocada na seguinte forma por meio da inversão das linhas e colunas; aqui, B e F são quadradas e 0 é uma matriz nula. 3ISSAI SCHUR (1875–1941), matemático alemão, também conhecido por seus importantes trabalhos em teoria de grupos. 4OSKAR PERRON (1880–1975), GEORG FROBENIUS (1849–1917), LOTHAR COLLATZ (1910–1990), matemáticos alemães, conhecidos por seus trabalhos em teoria do potencial, EDOs (Seção 5.4) e teoria de grupos e teoria numérica, respectivamente. 76 Parte E • Análise Numérica B C0 F O teorema de Perron tem várias aplicações, por exemplo, em economia. É interessante que ele pode ser obtido a partir de um teorema que fornece um algoritmo numérico: TEOREMA 6 Teorema da Inclusão de Collatz4 Consideremos que A = [ajk] seja uma matriz real n n cujos elementos são todos positivos. Chamemos de x um vetor real cujas componentes x1, • • • , xn são positivas e chamemos de y1, • • • , yn as componentes do vetor y = Ax. Então, o intervalo fechado no eixo real limitado pelo menor e pelo maior dos n coeficientes qj = yj/xj contém pelo menos um autovalor de A. PROVA Temos Ax = y ou (6) y Ax 0. A transposta AT satisfaz às condições do Teorema 5. Logo, AT possui um autovalor positivo l e, correspondente a esse autovalor, um autovetor u cujas componentes uj são todas positivas. Portanto, ATu = lu, e transpondo, obtemos uTA = luT. Disso e de (6), temos uT(y Ax) uTy uTAx uTy luTx uT(y lx) 0 ou escrevendo de forma compacta n j1 uj(yj lxj) 0. Visto que todas as componentes uj são positivas, segue-se que (7) yj lxj 0, yj lxj 0, ou seja, ou seja, qj l qj l para ao menos um j, para ao menos um j. e Como A e AT têm os mesmos autovalores, l é um autovalor de A e, de (7), decorre o enunciado do teorema. EXEMPLO 4 Limites para a Obtenção de Autovalores pelo Teorema de Collatz. Iteração Para uma dada matriz A com elementos positivos, escolhemos um x = x0 e efetuamos iterações, ou seja, computamos x1 = Ax0, x2 = Ax1, • • • , x20 = Ax19. Em cada passo, fazendo x = xj e y = Axj = xj + 1 calculamos um intervalo de inclusão pelo Teorema de Collatz. Isso fornece (6S) A 0,49 0,02 0,22 0,02 0,28 0,20 0,22 0,20 0,40 , x0 1 1 1 , x1 0,73 0,50 0,82 , x2 0,5481 0,3186 0,5886 , • • • , x19 0,00216309 0,00108155 0,00216309 , x20 0,00155743 0,000778713 0,00155743 e os intervalos 0,5 l 0,82, 0,3186/0,5 = 0,6372 l 0,5481/0,73 = 0,7508522 etc. Esses intervalos têm os comprimentos j 1 2 3 10 15 20 Tamanho 0,32 0,113622 0,0539835 0,0004217 0,0000132 0,0000004 Usando o polinômio característico, você pode verificar que os autovalores de A são 0,72, 0,36, 0,09, de modo que esses intervalos incluem o maior autovalor, 0,72. Seus comprimentos diminuíram com j, de modo que valeu a iteração foi proveitosa. A razão disso aparecerá na próxima seção, onde discutiremos um método iterativo para autovalores. P R O B L E M A S P R O P O S T O S 2 0 . 7 1–6 DISCOS DE GERSCHGORIN Encontre e esboce os discos ou intervalos que contêm os autovalores. Se você possui um aplicativo computacional à disposição, encontre o espectro e compare. 1. 6 1 0 1 7 1 0 1 8 Capítulo 20: Métodos Numéricos de Álgebra Linear 79 Se desejarmos uma seqüência convergente de autovetores, então no início de cada passo escalonamos o vetor, ou seja, dividimos suas componentes por aquela de maior valor absoluto, como no Exemplo 1 seguinte. EXEMPLO 1 Aplicação do Teorema 1. Ajuste de Escala No Exemplo 4 da Seção 20.7, para a matriz simétrica A e x0 = [1 1 1]T, obtemos de (1) e (2) e do escalonamento indicado A 0,49 0,02 0,22 0,02 0,28 0,20 0,22 0,20 0,40 , x0 1 1 1 , x1 0,890244 0,609756 1 , x2 0,931193 0,541284 1 x5 0,990663 0,504682 1 , x10 0,999707 0,500146 1 , x15 0,999991 0,500005 1 . Aqui, Ax0 = [0,73 0,5 0,82]T, com ajuste de escala para x1 = [0,73/0,82 0,5/0,82 1]T etc. O autovalor dominante é 0,72, e um autovetor é [1 0,5 1]T. Os q e δ correspondentes são computados antes de cada escalonamento. Portanto, no primeiro passo, q m1 m0 x0 TAx0 x0 Tx0 2,05 3 0,683333 d m2m0 q2 1/2 (Ax0)TAx0x0Tx0 q2 1/2 1,45533 q21/2 0,134743. Isso fornece os seguintes valores de q, δ e do erro e = 0,72 – q (calculados com 10D e arredondados para 6D): j 1 2 5 10 q 0,683333 0,716048 0,719944 0,720000 d 0,134743 0,038887 0,004499 0,000141 e 0,036667 0,003952 0,000056 5 108 Os limites dos erros são muito maiores que os erros reais. Isso é usual, embora os limites não possam ser melhorados; isto é, para matrizes simétricas especiais, os limites equivalem-se aos erros. Nossos atuais resultados são um pouco melhores que os obtidos pelo método de Collatz no Exemplo 4 da Seção 20.7, porém à custa de um número maior de operações. O deslocamento espectral, ou seja, a transição de A para A – kI, desloca todo autovalor por um fator –k. Embora achar um processo automático de obtenção de um bom k seja algo difícil, pode ser de alguma valia a utilização de algum outro método ou de pequenos experimentos computacionais preliminares. No Exemplo 1, o teorema de Gerschgorin fornece –0,02 λ 0,82 para o espectro inteiro (verifique!). O deslocamento por um fator –0,4 pode ser excessivo (pois então –0,42 λ 0,42), de modo que tentamos –0,2. EXEMPLO 2 Método de Potências com Deslocamento Espectral Para A – 0,2I com A como no Exemplo 1, obtemos as seguintes melhorias substanciais (onde o índice 1 se refere ao Exemplo 1 e o índice 2 ao exemplo atual). j 1 2 5 10 d1 0,134743 0,038887 0,004499 0,000141 d2 0,134743 0,034474 0,000693 1,8 106 e1 0,036667 0,003952 0,000056 5 108 e2 0,036667 0,002477 1,3 106 9 1012 P R O B L E M A S P R O P O S T O S 2 0 . 8 1–7 MÉTODO DE POTÊNCIAS COM ESCALONAMENTO Aplique o método de potências (3 passos) com escalonamento, usando x0 = [1 1]T ou [1 1 1]T ou [1 1 1 1]T, conforme a aplicabilidade. Forneça os coeficientes de Rayleigh e os limites dos erros. Mostre os detalhes do que fizer. 80 Parte E • Análise Numérica 20.9 Tridiagonalização e Fatoração QR Consideremos o problema de calcular todos os autovalores de uma matriz real simétrica A = [ajk], discutindo um método extensamente usado na prática. No primeiro estágio, reduzimos a matriz dada a uma matriz tridiagonal, ou seja, a uma matriz com todos os elementos diferentes de zero na diagonal principal e nas posições imediata- mente adjacentes à diagonal principal (como A3 na Fig. 447, Terceiro Passo). Essa redução foi concebida por A. S. Householder (J. Assn. Comput. Machinery 5 (1958), 335–342). Veja também a Ref. [E29] no Apêndice 1. Essa tridiagonalização de Householder simplificará a matriz sem alterar seus autovalores. Os últimos serão então determinados (aproximadamente) fatorando-se a matriz tridiagonalizada, conforme discutiremos depois nesta seção. Método da Tridiagonalização de Householder Dada uma matriz n n real simétrica A = [ajk], reduzimo-la por meio de n – 2 sucessivas transformações de similaridade (veja a Seção 20.6) envolvendo as matrizes P1, • • • , Pn–2 à forma tridiagonal. Essas matrizes são ortogonais e simétricas. Portanto, P1–1 = P1T = P1 e algo similar ocorre com as outras. Essas transformações produ- zem, a partir de uma dada A0 = A = [ajk], as matrizes A1 = [ajk(1)], A2 = [ajk(2)], • • • , An-2 = [ajk(n–2)] na forma A1 P1A 0P1 A2 P2A1P2 (1) • • • • • • • • • • • • B An2 Pn2An3Pn2. As transformações (1) criam os zeros necessários, no primeiro passo na Linha 1 e Coluna 1, no segundo passo na Linha 2 e Coluna 2 etc., como ilustra a Fig. 447 para a matriz 5 5. B é tridiagonal. Como podemos determinar P1, P2, • • • , Pn–2? Agora, todos esses Pr são da forma (2) Pr I 2vrvrT (r 1, • • • , n 2) 1. 3,5 2,02,0 0,5 2. 0,6 0,8 0,8 0,6 3. 5 22 7 4. 2 2 3 2 1 6 3 6 2 5. 2 1 1 1 3 2 1 2 3 6. 0 4 0 1 4 1 2 8 0 2 3 2 1 8 2 2 7. 5 1 0 0 1 3 1 0 0 1 3 1 0 0 1 5 8. (δ ótimo) No Problema 2, escolha x0 = [3 –1]T e mostre que q = 0 e δ = 1 para todos os passos, e que os autovalores são ±1, de modo que o intervalo [q – δ, q + δ] não pode em geral ser encurtado! Experimente com outros x0. 9. Prove que, em (2), se x é um autovetor, então δ = 0. Dê dois exemplos. 10. (Quociente de Rayleigh) Por que q é em geral uma aproximação do autovalor de maior valor absoluto? Quando q será uma boa aproximação? 11. (Deslocamento espectral, menor autovalor) No Problema 5, faça B = A – 3I (como talvez sugerido pelos elementos da diagonal) e verifique se você pode obter uma seqüência de q’s convergindo para um autovalor de A que seja o menor deles (não o maior) em valor absoluto. Use x0 = [1 1 1]T. Execute 8 passos. Verifique que A tem o espectro {0, 3, 5}. 12. EXPERIMENTO DE ÁLGEBRA COMPUTACIONAL. Méto- do de Potências com Escalonamento. Deslocamento. (a) Escre- va um programa para matrizes n n que imprima cada passo. Aplique-o à matriz (não-simétrica!) (20 passos), começando de [1 1 1]T. A 15 18 19 12 44 36 3 18 7 . (b) Experimente em (a) com o deslocamento. Qual é o desloca- mento ótimo encontrado? (c) Escreva um programa como em (a), porém para matrizes simé- tricas, e que imprima os vetores, os vetores escalonados, q, e δ. Aplique-o à matriz do Problema 6. (d) Encontre uma matriz (não-simétrica) para a qual δ em (2) já não seja mais um limite do erro. (e) Faça experimentos sistemáticos com a velocidade de conver- gência, escolhendo matrizes com o segundo maior autovalor (i) quase igual ao maior, (ii) um pouco diferente e (iii) muito dife- rente. Capítulo 20: Métodos Numéricos de Álgebra Linear 81 onde I é a matriz-identidade n n e vr = [vjr] é um vetor unitário com as suas r primeiras componentes iguais a 0; portanto, (3) v1 0 • • • , v2 0 0 • • • , • • • , vn2 0 0 • • • onde os asteriscos denotam as outras componentes (que serão diferentes de zero em geral). Passo 1. v1 tem as componentes (4) v11 0 (a) v21 1 (b) vj1 onde (c) S1 a21 2 a31 2 • • • an1 2 aj1 sgn a21 2v21S1 a21 S1 1 2 j 3, 4, • • • , n onde S1 > 0, e sgn a21 = +1 se a21 0 e sgn a21 = –1 se a21 < 0. Com isso, calculamos P1 por (2) e então A1 por (1). Esse foi o primeiro passo. * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * Primeira Etapa Segunda Etapa Terceira Etapa A 1 = P 1 AP 1 A 2 = P2 A1P2 A3 = P3A2P3 Fig. 447. Método de Householder para uma matriz 5 5. Os espaços em branco são zeros criados pelo método Passo 2. Calculamos v2 por (4) com todos os subscritos aumentados de 1 unidade e os ajk substituídos por a jk(1), que são os elementos de A1 recém-calculados. Assim [veja também (3)], (4*) v12 v22 0 v32 1 vj2 S2 a32 (1)2 a42 (1)2 • • • an2 (1)2 . a j2 (1) sgn a32 (1) 2v32S2 a32 (1) S2 1 2 j 4, 5, • • • , n onde Com isso, calculamos P2 por (2) e então A2 por (1). Passo 3. Calculamos v3 por (4*) com todos os subscritos aumentados de 1 unidade e os a jk(1) substituídos pelos elementos a jk(2) de A2, e assim por diante. 84 Parte E • Análise Numérica Como Obter a Fatoração QR, ou seja, B = B0 = [bjk] = Q0R0. A matriz tridiagonal B tem geralmente n – 1 elementos diferentes de zero abaixo da diagonal principal. Eles são b21, b32, • • • , bn,n–1. A partir da esquerda, multiplicamos B por uma matriz C2 tal que C2B = [b jk(2)] tenha b21(2) = 0. Multiplicamos isso por uma matriz C3 tal que C3C2B = [b jk(3)] tenha b32(3) = 0 etc. Após n – 1 multiplicações dessas, chegamos a uma matriz triangular superior R0, a saber, (7) CnCn1 • • • C3C2B0 R0. Essas matrizes n n Cj são muito simples. Cj tem a submatriz 2 2 cos ujsen uj sen uj cos uj (uj apropriado) em Linhas j – 1 e j e Colunas j – 1 e j; nos demais lugares da diagonal principal, a matriz Cj tem elementos 1; e todos os outros elementos são 0. (Essa submatriz é a matriz de uma rotação de uma rotação plana por um ângulo θj ; veja o Projeto de Equipe 28 da Seção 7.2.) Por exemplo, se n = 4, escrevendo cj = cos θj e sj = sen θj, temos C2 c2 s2 0 0 s2 c2 0 0 0 0 0 0 1 0 0 1 , C3 1 0 0 0 0 c3 s3 0 0 s3 c3 0 0 0 0 1 , C4 1 0 0 1 0 0 0 0 0 0 c4 s4 0 0 s4 c4 . Essas Cj são ortogonais. Logo, seu produto em (7) é ortogonal e o mesmo ocorre com a inversa desse produto. Chamemos essa inversa de Q0. Então, de (7), (8) B0 Q0R0 onde, com Cj–1 = CjT, (9) Q0 (CnCn1 • • • C3C2)1 C2TC3T • • • Cn1TCnT. Essa é a nossa fatoração QR de B0. Disso, temos por (5) com s = 0 (10) B1 R0Q0 R0C2TC3T • • • Cn1TCnT. Não precisamos explicitamente de Q0, mas, para obter B1 de (10), primeiramente calculamos R0C2T, então (R0C2T) C3T etc. Algo similar ocorre nos demais passos, que produzem B2, B3, • • •. Determinação de cos θj e sen θj. Finalmente, mostremos como encontrar os ângulos de rotação. Em C2, cos θ2 e sen θ2 devem ser tais que b21(2) = 0 no produto C2B c2 s2 • • s2 c2 • • 0 0 • • • • • • • • • • • • • • b11 b21 • • b12 b22 • • b13 b23 • • • • • • • • • • • • • • . Agora, b21(2) é obtida multiplicando-se a segunda linha de C2 pela primeira coluna de B, b21(2) s2b11 c2b21 (sen u2)b11 (cos u2)b21 0. Logo, tan θ2 = s2/c2 = b21/b11, e (11) cos 2 sen 2 . b21/b11 1 (b21/b11) 2 tan 2 1 tan2 2 1 1 (b21/b11) 2 1 1 tan2 2 Algo similar ocorre para θ3, θ4, • • • O próximo exemplo ilustra tudo isso. Capítulo 20: Métodos Numéricos de Álgebra Linear 85 EXEMPLO 2 Método da Fatoração QR Calcule todos os autovalores da matriz A 6 4 1 1 4 6 1 1 1 1 5 2 1 1 2 5 . Solução. Primeiramente, reduzimos A à forma tridiagonal. Aplicando o método de Householder, obtemos (veja o Exemplo 1) A2 6 18 0 0 18 7 2 0 0 2 6 0 0 0 0 3 . Do determinante característico vemos que A2, e conseqüentemente A, tem o autovalor 3. (Você pode ver isso diretamente de A2?) Logo, isso é suficiente para a aplicação do método QR à matriz tridiagonal 3 3 B0 B 6 18 0 18 7 2 0 2 6 . Passo 1. Multiplicamos B à esquerda por C2 cos u2 sen u2 0 sen u2 cos u2 0 0 0 1 e então, C2B por C3 1 0 0 0 cos u3 sen u3 0 sen u3 cos u3 . Aqui, (–sen θ2) 6 + (cos θ2)(– 18 ) = 0 fornece (11) cos θ2 = 0,816 496 58 e sen θ2 = –0,577 350 27. Com esses valores, calculamos C2B 7,348 469 23 0 0 7,505 553 50 3,265 986 32 1,414 213 56 0,816 496 58 1,154 700 54 6,000 000 00 . Em C3, temos de (–sen θ3) 3,265 986 32 + (cos θ3) 1,414 213 56 = 0 os valores cos θ3 = 0,917 662 94 e sen θ3 = 0,397 359 71. Isso fornece R0 C3C2B 7,348 469 23 0 0 7,505 553 50 3,559 026 08 0 0,816 496 58 3,443 784 13 5,047 146 15 . Disso, calculamos B1 R0C2 TC3 T 10,333 333 33 2,054 804 67 0 2,054 804 67 4,035 087 72 2,005 532 51 0 2,005 532 51 4,631 578 95 que é simétrica e tridiagonal. Os elementos fora da diagonal em B1 são ainda grandes em valor absoluto. Logo, temos ainda que continuar. Passo 2. Fazemos os mesmos cálculos que no primeiro passo, com B0 = B substituído por B1 e C2 e C3 adequadamente modificados, com os novos ângulos sendo θ2 = –0,196 291 533 e θ3 = 0,513 415 589. Obtemos então, R1 10,535 653 75 0 0 2,802 322 41 4,083 295 84 0 0,391 145 88 3,988 240 28 3,068 326 68 e, disso, B2 10,879 879 88 0,796 379 18 0 0,796 379 18 5,447 386 64 1,507 025 00 0 1,507 025 00 2,672 733 48 . 86 Parte E • Análise Numérica Vemos que os elementos de fora da diagonal são menores em valor absoluto que os de B1, porém são ainda muito grandes para permitir que os elementos da diagonal sejam boas aproximações dos autovalores de B. Passos Adicionais. Listamos os elementos da diagonal principal e o elemento de fora da diagonal com o maior valor absoluto, que é b12 ( j) b21 ( j) em todos os passos. Você pode mostrar que a matriz A dada possui o espectro 11, 6, 3, 2. Passo j b11 ( j) b22 ( j) b33 ( j) max jk b jk (J) 3 10,966 892 9 5,945 898 56 2,087 208 51 0,585 235 82 5 10,997 087 2 6,001 815 41 2,001 097 38 0,120 653 34 7 10,999 742 1 6,000 244 39 2,000 013 55 0,035 911 07 9 10,999 977 2 6,000 022 67 2,000 000 17 0,010 684 77 Vendo nossa discussão em retrospecto, reconhecemos que o propósito da aplicação da tridiagonalização de House- holder antes do método da fatoração QR é a redução substancial do custo em cada fatoração QR, em particular se A for grande. A aceleração da convergência e, portanto, a redução adicional do custo pode ser obtida por um deslocamento espectral, isto é, tomando-se Bs – ksI ao invés de Bs com um ks adequado. Possíveis escolhas para ks são discu- tidas na Ref. [E29]. P R O B L E M A S P R O P O S T O S 2 0 . 9 1–4 TRIDIAGONALIZAÇÃO DE HOUSEHOLDER Tridiagonalize, mostrando os detalhes: 1. 3,5 1,0 1,5 1,0 5,0 3,0 1,5 3,0 3,5 2. 0 1 1 1 0 1 1 1 0 3. 0,98 0,04 0,44 0,04 0,56 0,40 0,44 0,40 0,80 4. 8 8 2 2 8 8 2 2 2 2 6 4 2 2 4 6 5–9 FATORAÇÃO QR Execute três passos da fatoração QR para encontrar aproximações dos autovalores da: 5. Matriz na resposta do Problema 1 6. Matriz na resposta do Problema 3 7. 18 2 0 2 0 8 2 2 2 8. 16,2 0,1 0 0,1 4,3 0,2 0 0,2 4,1 9. 7,0 0,1 0 0,1 4,0 0,1 0 0,1 1,0 10. EXPERIMENTO DE ÁLGEBRA COMPUTACIONAL. Méto- do QR. Tente encontrar experimentalmente de quais propriedades da matriz a velocidade de decrescimento dos elementos de fora da diagonal depende. Para esse propósito, escreva um programa que primeiramente tridiagonalize e então execute os passos QR. Aplique o programa às matrizes dos problemas 1, 3, e 4. Sintetize seus resultados em um pequeno texto. Q U E S T Õ E S E P R O B L E M A S D E R E V I S Ã O D O C A P Í T U L O 2 0 1. Quais são as principais áreas problemáticas na álgebra linear numé- rica? 2. O que é pivotação? Quando e como podemos aplicá-la? 3. O que acontece se aplicarmos a eliminação de Gauss a sistemas que não possuem solução? 4. O que é o método de Doolittle? E qual a sua relação com a elimi- nação de Gauss? 5. O que é o método de Cholesky? Quando podemos aplicá-lo? 6. O que você sabe sobre a convergência do método de Gauss– Seidel? 7. O que é mau condicionamento? O que é o número condicional e o seu significado? 8. O que é a aproximação por mínimos quadrados? O que são as equações normais? 9. O que é um autovalor de uma matriz? Por que os problemas de autovalores são importantes? Dê alguns exemplos típicos. C A P Í T U L O 2 1 Métodos Numéricos para Equações Diferenciais Ordinárias (EDOs) e Equações Diferenciais Parciais (EDPs) Os métodos numéricos para resolver equações diferenciais são de grande importância prática em engenha- ria e física, onde os problemas práticos freqüentemente levam a equações diferenciais que não podem ser resolvidas pelos métodos vistos nos Capítulos. de 1 a 6, nem pelos vistos no Capítulo 12, nem por métodos similares. Além disso, às vezes ocorre que uma EDO chega, de fato, a ter uma fórmula de solução (como as EDOs das Seções 1.3–1.5); a qual, entretanto, em alguns casos específicos, pode se tornar tão complicada que é preferível aplicar um método numérico para se resolver a equação. Este capítulo explica e fornece aplicações dos métodos básicos para se obterem soluções numéricas de EDOs (Seções 21.1–21.3) e EDPs (Seções 21.4–21.7). As Seções 21.1 e 21.2 podem ser estudadas imediatamente após o Capítulo 1, e a Seção 21.3 imedia- tamente após o Capítulo 2, visto que estas seções independem dos Capítulos 19 e 20. As Seções 21.4–21.7, sobre EDPs, podem ser estudadas imediatamente após o Capítulo 12 se os estu- dantes já tiverem algum conhecimento de sistemas lineares de equações algébricas. Pré-requisito: Seções 1.1–1.5 para EDOs, Seções 12.1–12.3, 12.5, 12.10 para EDPs. Referências e Respostas dos Problemas: Parte E do Apêndice 1 (veja também as Partes A e C) e Apên- dice 2. 21.1 Métodos para EDOs de Primeira Ordem Do Capítulo 1, sabemos que uma EDO de primeira ordem tem a forma F(x, y, y) = 0 e freqüentemente pode ser escrita na forma explícita y = f(x, y). Um problema de valor inicial para essa equação tem a forma (1) y ƒ(x, y), y(x0) y0 onde x0 e y0 são dados e supomos que o problema tenha uma solução única em algum intervalo aberto a < x < b contendo x0. Nesta seção, discutiremos métodos para calcular valores numéricos aproximados da solução y(x) de (1) em pontos eqüidistantes sobre o eixo x x1 x0 h, x2 x0 2h, x3 x0 3h, • • • onde o tamanho do passo h é um número fixo, por exemplo, 0,2 ou 0,1 ou 0,001, cuja escolha discutiremos mais tarde nesta seção. Esses são os métodos passo a passo, que usam uma mesma fórmula em cada um dos passos. Essas fórmulas são sugeridas pela série de Taylor (2) y(x h) y(x) hy(x) h2 2 y(x) • • • . Para um h pequeno, as potências mais elevadas h2, h3, • • • dão resultados muito pequenos. Isto sugere a aproxi- mação grosseira y(x h) y(x) hy(x) y(x) hƒ(x, y) (com a segunda linha obtida da EDO dada) e o seguinte processo de iteração. No primeiro passo, calculamos 90 Parte E • Análise Numérica y1 = y0 + hf(x0, y0) que é uma aproximação de y(x1) = y(x0 + h). No segundo passo, computamos y2 = y1 + hf(x1, y1) que é uma aproximação de y(x2) = y(x0 + 2h) etc. e, em geral, (3) yn1 yn hƒ(xn, yn) (n 0, 1, • • •). Este é o chamado método de Euler ou de Euler–Cauchy. Geometricamente, ele é uma aproximação da curva y(x) por um polígono cujo primeiro lado é tangente a esta curva em x0 (veja a Fig. 448). y x hh Inclinação f(x 0 , y 0 ) Inclinação f(x 1 , y 1 ) y 0 y 1 y 2 x 0 x 1 x 2 Fig. 448. Método de Euler Este método grosseiro é raramente usado na prática mas, devido à sua simplicidade, ele explica bem o princípio dos métodos baseados na série de Taylor. A fórmula de Taylor com um resto se escreve como y(x h) y(x) hy(x) 1_2h2y(j) (onde x j x h). Ela mostra que, no método de Euler, o erro de truncamento em cada passo, ou erro de truncamento local, é proporcional a h2 e se escreve como O(h2), onde O sugere ordem (veja também a Seção 20.1). Ora, em um intervalo fixo em x onde desejamos resolver uma EDO, o número de passos é proporcional a 1/h. Logo, o erro total ou erro global é proporcional a h2(1/h) = h1. Por este motivo, o método de Euler é chamado de método de primeira ordem. Além disso, há erros de arredondamento neste e em outros métodos, que podem afetar a precisão dos valores y1, y2, • • • progressivamente à medida que n cresce, conforme veremos. Tabela 21.1 Método de Euler Aplicado a (4) no Exemplo 1 e o Erro n xn yn 0,2(xn yn) Valores exatos Erro en 0 0,0 0,000 0,000 0,000 0,000 1 0,2 0,000 0,040 0,021 0,021 2 0,4 0,040 0,088 0,092 0,052 3 0,6 0,128 0,146 0,222 0,094 4 0,8 0,274 0,215 0,426 0,152 5 1,0 0,489 0,718 0,229 EXEMPLO 1 Método de Euler Aplique o método de Euler ao seguinte problema de valor inicial, escolhendo h = 0,2 e calculando y1, • • • , y5 (4) y x y, y(0) 0. Solução. Aqui f(x, y) = x + y; logo f (xn, yn) = xn + yn, e vemos que (3) se torna yn1 yn 0,2(xn yn). A Tabela 21.1 mostra os cálculos e os valores da solução exata para y(x) = ex – x – 1 obtida a partir de (4) na Seção 1.5, e os erros respectivos. Na prática, a solução exata é desconhecida, mas pode-se ter uma idéia da precisão dos valores aplicando-se mais uma vez o método de Euler, com a iteração 2h = 0,4, fazendo yn* representar a aproximação agora obtida e comparando-se as aproximações correspondentes. Este cálculo é: Capítulo 21: Métodos Numéricos para Equações Diferenciais Ordinárias (EDOs) e Equações Diferenciais Parciais (EDPs) 91 xn yn* 0,4(xn yn) yn na Tabela 21.1 Diferença yn yn* 0,0 0,000 0,000 0,000 0,000 0,4 0,000 0,160 0,040 0,040 0,8 0,160 0,274 0,114 Chamemos de en e en* os erros dos cálculos com h e 2h, respectivamente. Uma vez que o erro é de ordem h2, quando trocamos h por 2h, esse erro fica multiplicado por 22 = 4, mas, visto que necessitamos somente de metade dos passos de antes, ele será multiplicado somente por 4/2 = 2. Logo, en* 2en de modo que a diferença é en* en 2en en en. Ora, y yn en yn* en* pela definição de erro; por conseguinte, en* en yn yn* indica en qualitativamente. Em nossos cálculos, y2 – y2* = 0,04 – 0 = 0,04 (o erro real é 0,052, veja a Tabela 21.1) e y4 – y4* = 0,274 – 0,160 = 0,114 (na verdade, 0,152). EXEMPLO 2 Método de Euler para uma EDO Não-linear A Fig. 449 diz respeito ao problema de valor inicial (5) y (y 0,01x2)2 sen (x2) 0,02 x, y(0) 0,4 e mostra a curva da solução y = 1/[2,5 – S(x)] + 0,01x2, onde S(x) é a integral de Fresnel (38) no Apêndice 3.1. Ela também mostra 80 valores aproximados para 0 x 4 obtidos pelo método de Euler a partir de (3), yn1 yn 0,05 [(yn 0,01xn2)2 sen (xn2) 0,02 xn] . Embora h = 0,05 seja menor do que h no Exemplo 1, a precisão ainda não é boa. É interessante notar que o erro não aumenta monotonamente, o que é óbvio, visto que a solução não é monótona. Retornaremos a esta EDO em Problemas Propostos. y x0 2 41 3 0,40 0,50 0,60 0,70 Fig. 449. Curva de solução e aproximação de Euler do Exemplo 2 Seleção Automática do Tamanho do Passo Variável Usando-se Modernos Programas Computacionais A idéia de integração adaptativa, que foi justificada e explicada na Seção 19.5, se aplica igualmente bem à solução de EDOs por métodos numéricos. Ela é agora utilizada para alterar automaticamente o tamanho h do passo, dependendo da variabilidade de y = f determinado por (6*) y ƒ ƒx ƒyy ƒx ƒyƒ. Dessa forma, um programa computacional moderno seleciona automaticamente os tamanhos hn dos passos de tal modo que o erro da solução não exceda um tamanho máximo dado TOL (termo que sugere tolerância). Ora, para o método de Euler, quando o tamanho do passo é h = hn, o erro local em xn vale cerca de 12hn 2 y(xn). Nossa exigência é de que ele seja igual a uma dada tolerância TOL, (6) (a) 12hn 2y(xn) TOL, portanto (b) hn 2 TOL y(xn) . y(x) não pode ser nula no intervalo J: x0 x = xN no qual se deseja obter uma solução. Chamemos de K o mínimo valor de y(x) em J e suponhamos que K > 0. Por (6), o mínimo valor de y(x) corresponde ao máximo valor de h H 2 TOL/K de (6). Portanto, 2 TOL HK. Podemos inserir isto em (6b), obtendo diretamente por álgebra, (7) hn w(xn)H onde w(xn) Ky(xn) . 94 Parte E • Análise Numérica aqui, f é tal que este problema tem uma solução única no intervalo [x0, xN] (veja a Seção 1.7). ENTRADA: Função f, valores iniciais x0, y0, tamanho do passo h, número de iterações N SAÍDA: Aproximação yn+1 da solução y(xn+1) em xn+1 = x0 + (n + 1)h, onde n = 0, • • • , N – 1 Para n 0, 1, • • • , N 1 faça: k1 hƒ(xn, yn) k2 hƒ(xn 1_ 2h, yn 1_ 2k1) k3 hƒ(xn 1_ 2h, yn 1_ 2k2) k4 hƒ(xn h, yn k3) xn 1 xn h yn 1 yn 1_ 6(k1 2k2 2k3 k4) SAÍDA xn 1, yn 1 Fim Pare Fim RUNGE–KUTTA EXEMPLO 4 Método Clássico de Runge–Kutta Aplique o método de Runge–Kutta ao problema de valor inicial(4) no Exemplo 1, escolhendo h = 0,2 como antes e fazendo cinco passos. Solução. Para este problema, temos f(x, y) = x + y. Logo, k1 0,2(xn yn), k2 0,2(xn 0,1 yn 0,5k1), k3 0,2(xn 0,1 yn 0,5k2), k4 0,2(xn 0,2 yn k3). A Tabela 21.5 mostra os resultados e seus erros, que são menores por fatores de 103 e 104 do que aqueles para os dois métodos de Euler. Veja também a Tabela 21.6. Mencionamos de passagem que, como os presentes k1, k2, k3, k4 são simples, poderíamos poupar operações substituindo k1 por k2 e então k2 por k3 etc.; a fórmula resultante é mostrada na Coluna 4 da Tabela 21.5. Tabela 21.5 Método de Runge–Kutta Aplicado a (4) n xn yn 0,2214(xn yn) 0,0214 Valores exatos (6D) y ex x 1 106 Erro de yn 0 0,0 0 0,021 400 0,000 000 0 1 0,2 0,021 400 0,070 418 0,021 403 3 2 0,4 0,091 818 0,130 289 0,091 825 7 3 0,6 0,222 107 0,203 414 0,222 119 12 4 0,8 0,425 521 0,292 730 0,425 541 20 5 1,0 0,718 251 0,718 282 31 Tabela 21.6 Comparação da Precisão dos Três Métodos Considerados no Caso do Problema de Valor Inicial (4), com h = 0,2 Erro x y ex x 1 Euler (Tabela 21.1) Euler aperfeiçoado (Tabela 21.3) Runge–Kutta (Tabela 21.5) 0,2 0,021 403 0,021 0,0014 0,000 003 0,4 0,091 825 0,052 0,0034 0,000 007 0,6 0,222 119 0,094 0,0063 0,000 011 0,8 0,425 541 0,152 0,0102 0,000 020 1,0 0,718 282 0,229 0,0156 0,000 031 Capítulo 21: Métodos Numéricos para Equações Diferenciais Ordinárias (EDOs) e Equações Diferenciais Parciais (EDPs) 95 Erro e Controle do Tamanho do Passo. RKF (Runge–Kutta–Fehlberg) A idéia da integração adaptativa (Seção 19.5) tem casos análogos para os métodos de Runge–Kutta (e outros). Na Tabela 21.4 sobre RK (Runge–Kutta), se em cada passo calcularmos as aproximações y e y com tamanhos do passo h e 2h, respectivamente, ocorre um erro por iteração igual a 25 = 32 vezes o do anterior; todavia, como temos apenas a metade do número de passos para 2h, o fator real é 25/2 = 16, de modo que, por exemplo, e(2h) 16e(h) e portanto y(h) y(2h) e(2h) e(h) (16 1)e(h). Logo, o erro e e(h) para o tamanho do passo é cerca de (10) e 115 (y y ) onde y y y(h) y(2h), como se disse antes. A Tabela 21.7 ilustra (10) para o problema de valor inicial (11) y (y x 1)2 2, y(0) 1, o tamanho do passo h = 0,1 e 0 x 0,4. Vemos que a estimativa está próxima do erro real. Esse método de estimativa de erro é simples, mas pode ser instável. Tabela 21.7 Método de Runge–Kutta Aplicado ao Problema de Valor Inicial (11) e Estimativa de Erro (10). Solução Exata y = tan x + x + 1 x y (Tamanho da iteração h) y (Tamanho da iteração 2h) Erro estimado (10) Erro real Solução exata (9D) 0,0 1,000 000 000 1,000 000 000 0,000 000 000 0,000 000 000 1,000 000 000 0,1 1,200 334 589 0,000 000 083 1,200 334 672 0,2 1,402 709 878 1,402 707 408 0,000 000 165 0,000 000 157 1,402 710 036 0,3 1,609 336 039 0,000 000 210 1,609 336 250 0,4 1,822 792 993 1,822 788 993 0,000 000 267 0,000 000 226 1,822 793 219 RKF. E. Fehlberg [Computing 6 (1970), 61–71] propôs e desenvolveu o controle do erro pelo uso de dois méto- dos RK de diferentes ordens para ir de (xn, yn) a (xn+1, yn+1). A diferença entre os valores calculados de y em xn+1 fornece uma estimativa de erro a ser usada para o controle do tamanho do passo. Fehlberg descobriu duas fórmulas RK que, juntas, precisam apenas de seis resultados por iteração. Apresentamos estas fórmulas aqui porque o RKF tornou-se bem popular. Por exemplo, o Maple o utiliza (também para sistemas de EDOs). O método RK de quinta ordem de Fehlberg é (12a) yn1 yn g1k1 • • • g6k6 com o vetor coeficiente g [g1 • • • g6], (12b) g [ 16135 0 665612825 2856156430 950 255 ]. Seu método RK de quarta ordem é (13a) y*n1 yn g1*k1 • • • g5*k5 com o vetor coeficiente (13b) g* [ 25216 0 14082565 21974104 15]. Em ambas as fórmulas, usamos simultaneamente por iteração somente seis resultados, quais sejam, (14) _11 40k5). _845 4104k4) _1859 4104k4 _7296 2197k3) _3680 513 k3 _3544 2565k3 _9 32k2) _7200 2197k2) 8k2) 2k2) 1_ 4k1) _3 32k1) _1932 2197k1) _439 216k1) _8 27k1) yn yn yn yn yn k1 hƒ(xn, yn) k2 hƒ(xn 1_ 4h, k3 hƒ(xn 3_ 8h, k4 hƒ(xn _12 13h, k5 hƒ(xn h, k6 hƒ(xn 1_ 2h, 96 Parte E • Análise Numérica A diferença entre (12) e (13) dá o erro estimado (15) n 1 ≈ yn 1 y*n 1 _1360k1 _1284275k3 _219775240k4 _150k5 _255k6. EXEMPLO 5 Runge–Kutta–Fehlberg Para o problema de valor inicial (11), obtemos de (12)–(14) com h = 0,1 na primeira iteração, os valores 12S k1 0,200000 000000 k3 0,200140 756867 k5 0,201006 676700 k2 0,200062 500000 k4 0,200856 926154 k6 0,200250 418651 y1* 1,200334 66949 y1 1,200334 67253 e o erro estimado e1 ≈ y1 y*1 0,000000 00304. O valor exato 12S é y(0,1) = 1,200334 67209. Logo, o erro real de y1 é –4,4 10–10, ou seja, menor por um fator de 200 que o da Tabela 21.7. A Tabela 21.8 resume as características essenciais dos métodos nesta seção. Pode-se mostrar que esses métodos são numericamente estáveis (definição na Seção 19.1). Eles são métodos de uma iteração, porque, em cada passo, utilizamos os dados de apenas uma iteração precedente, em contraste com os métodos multi-iterativos, onde em cada iteração utilizam-se dados de várias iterações anteriores, como veremos na próxima seção. Tabela 21.8 Métodos Considerados e Suas Ordens (= Seu Erro Global) Método Avaliação da função por iteração Erro global Erro local Euler 1 O(h) O(h2) Euler aperfeiçoado 2 O(h2) O(h3) RK (4a ordem) 4 O(h4) O(h5) RKF 6 O(h5) O(h6) Método de Euler Reverso. EDOs Rígidas A fórmula reversa de Euler para se obter uma solução numérica de (1) é (16) yn1 yn hƒ(xn1, yn1) (n 0, 1, • • •). Esta fórmula é obtida calculando-se o valor do lado direito da equação no novo ponto (xn+1, yn+1); este é o chamado esquema de Euler reverso. Para um valor conhecido de yn, ele fornece yn+1 implicitamente, de modo que define assim um método implícito, diferentemente do método de Euler (3), que fornece yn+1 explicitamente. Logo, é preciso resolver (16) para yn+1. O grau de dificuldade desse método depende de f em (1). Para uma EDO linear, não há qualquer problema, como mostra o Exemplo 6 (a seguir). Este método é particularmente útil para EDOs “rígidas”, como as que freqüentemente ocorrem no estudo de vibrações, circuitos elétricos, reações químicas, etc. A situação de rigidez é, grosso modo, a seguinte; detalhes podem, por exemplo, ser vistos em [E5], [E25] e [E26] no Apêndice 1. Nos métodos até aqui considerados, os termos referentes aos erros envolvem uma derivada de ordem superior. Poderíamos, então, nos perguntar o que acontece se fizermos h aumentar. Ora, se o erro (a derivada) aumentar rapidamente e a solução desejada também aumentar rapidamente, nada acontecerá. Entretanto, se essa solução não aumentar rapidamente, então, com o crescimento de h, o termo do erro pode atingir um tamanho tal que deixará o resultado completamente sem sentido, conforme se vê na Fig. 450. Portanto, chamamos de rígidos tanto a EDO para a qual h deve se restringir a pequenos valores, quanto o sistema físico modelado por essa EDO. O uso desse termo é sugerido pelo sistema massa-mola apresentando uma mola rígida (ou seja, uma mola com um grande valor de k; veja a Seção 2.4). O Exemplo 6 ilustra o fato de que os métodos implícitos removem a dificuldade decorrente do aumento de h no caso de rigidez; pode-se demonstrar que, na aplicação de um método implícito, a solução permanece estável sob qualquer aumento de h, embora a precisão diminua à medida que h aumenta.