Baixe o app para aproveitar ainda mais
Prévia do material em texto
INSTITUTO FEDERAL DE EDUCAÇÃO TECNOLÓGICA Av. Jerônimo Figueira da Costa, 3014 - Bairro Pozzobon. CEP 15503-110 - Votuporanga - SP www.ifsp.edu.br ENGENHARIA CIVIL: NOTAS DE AULA CÁLCULO NUMÉRICO Prof. Msc. Gustavo Cabrelli Nirschl Janeiro de 2018 Agradecimentos: Carlos Eduardo Alves da Silva - Técnico de TI do Campus Fernando Parreira - Técnico de TI do Campus Isabela Cassia Dominical Parra – orientada de iniciação científica relacionada aos assuntos desta apostila (2015 e 2016) e discente da 1ª turma do curso de Engenharia Civil Marcelo Yugi Nomoto – ajudou nos meus primeiros passos na programação web O QUE IMPORTA É O CONHECIMENTO E A COMUNICAÇÃO Caro estudante, o professor tem o dever não só de te passar conhecimento, mas de fazer você aprender a aprender. Isso é o que um bom professor deve ter como objetivo. Para isso, o professor te direciona aos estudos e avalia a tua evolução por meio das notas. E elas podem ser cruéis: são um número que não percebe se você está com problemas familiares, pessoais, de saúde. Naquele dia de avaliação você tem que ter aquele conteúdo na cabeça, ou o número será baixo! É o sistema! Você tem que se adequar! Como um futuro profissional, saiba avaliar a tua própria nota, individualmente, considerando o que está acontecendo com todo o resto da tua vida naquele momento. Se não for por preguiça de estudar, não se martirize e nem culpe o professor. Apenas refaça os estudos, repita o processo. Isso é a repetência! Isso não é vergonha! Quando você entrar no mercado de trabalho, o contratante irá te entrevistar para tentar perceber se você o ajudará a lucrar e a desenvolver a empresa, mesmo que não de imediato, mas num futuro próximo. Então ele irá tentar ver o que você sabe fazer, se aquilo é suficiente para tocar os trabalhos da empresa. Ele pode até querer ver tuas notas, mas elas não serão mais importantes do que o teu conhecimento adquirido. Hoje em dia ninguém consegue trabalhar isoladamente, sozinho num canto. Os grandes projetos exigem trocas de conhecimento entre diversas áreas, cooperação. É a multidisciplinaridade! Então as corporações irão gostar do teu trabalho se você souber se comunicar oralmente e de forma escrita. Se falar inglês, será um profissional diferenciado! Porque é bonito saber inglês? Não, porque os grandes projetos exigirão informações e estudos em inglês, terão contato com pessoas internacionais, e até serão executados em outros países. Saber escrever é mera formalidade? Também não! Você precisará de muitas pessoas para executar aquilo que está na tua cabeça e se você conseguir dizer a elas exatamente como deve ser feito, o projeto ficará bom. Não adianta saber tudo e não ter a chave da porta que leva tudo isso pro mundo exterior: a comunicação oral e escrita! Então, meu amigo e aluno, estude sim a técnica, e muito, mas estude português e inglês, no mínimo. Comunique-se com os colegas, não por amizade somente, mas para dar e receber novas informações a cada dia. Isso sim te levará a ser um ótimo Engenheiro Civil! Bons estudos! Prof. Gustavo Cabrelli Nirschl PROVAS E TRABALHOS DA DISCIPLINA 1º BIMESTRE: Media1 = (Prova1 * 0,8) + (Trabalhos1 * 0,2) 2º BIMESTRE: Media2 = (Prova2 * 0,8) + (Trabalhos2 * 0,2) MÉDIA = 0,4*Media1 + 0,6*Media2 SE 4,0 MÉDIA < 6,0 RECUPERAÇÃO FINAL CRITÉRIOS DE APROVAÇÃO FREQUENCIA ≥ 75% E MÉDIA FINAL ≥ 6,0 SOBRE ESTAS NOTAS DE AULA Estas notas de aula foram escritas, a partir de 2015, com base em outras bibliografias e na experiência acumulada do professor, para a disciplina de Cálculo Numérico do curso de Engenharia Civil do IFSP campus Votuporanga. Elas são constantemente revisadas devido a possíveis erros e necessidades dos alunos não identificadas previamente. Por isso é importante o aluno estar com a revisão mais recente e não reutilizar materiais de turmas anteriores. Lembra-se ao aluno que essas notas de aula não substituem os livros em hipótese alguma e não podem ser usadas para outro fim que não para os estudos desta disciplina. BIBLIOGRAFIA COMPLEMENTAR BARROSO, L.C. et al. Cálculo numérico: com aplicações. 2.ed. São Paulo: Harbra, 1987. LIPSCHUTZ, S. Álgebra Linear. 3.ed. São Paulo: Pearson, 1994. FARRER, H. et al. Pascal estruturado. 3.ed. Rio de Janeiro: LTC, 1999. FORBELLONE, A.L.; EBERSPACHER, H.F. Lógica de programação. 3.ed. São Paulo: Prentice Hall, 2005. RUGGIERO, M.A.G; LOPES, V.L.R. Cálculo numérico: aspectos teóricos e computacionais. 2.ed. São Paulo: Makron Books, 1996. SALIBA, L.W. Técnicas de programação: uma abordagem estruturada. São Paulo: Makron Books,1993. GUIDORIZZI, H.L. Um curso de cálculo. v.2. 5.ed.Rio de Janeiro: LTC, 2001. GUIDORIZZI, H.L. Um curso de cálculo. v.3. 5.ed.Rio de Janeiro: LTC, 2001. SUMÁRIO 1 MÉTODOS NUMÉRICOS NA SOLUÇÃO DE PROBLEMAS .................................. 6 2 INTRODUÇÃO À TEORIA DE ERRO E ESTABILIDADE ....................................... 9 3 SISTEMAS DE EQUAÇÕES LINEARES............................................................... 14 4 OPERAÇÕES COM MATRIZES ............................................................................ 40 5 ZEROS DE FUNÇÕES........................................................................................... 47 6 FUNÇÕES DE INTERPOLAÇÃO DE DADOS ...................................................... 57 7 FUNÇÕES DE APROXIMAÇÃO DE DADOS ........................................................ 61 8 INTEGRAÇÃO NUMÉRICA DE FUNÇÕES ........................................................... 72 9 DERIVAÇÃO NUMÉRICA DE FUNÇÕES ............................................................. 78 10 SOLUÇÃO DE EQUAÇÕES DIFERENCIAIS ...................................................... 83 11 CONSIDERAÇÕES FINAIS ................................................................................. 98 12 REFERÊNCIAS BIBLIOGRÁFICAS .................................................................. 103 APÊNDICE A - INTRODUÇÃO À PROGRAMAÇÃO EM LAZARUS/PASCAL ..... 105 APÊNDICE B - INTRODUÇÃO À PROGRAMAÇÃO EM HTML/JAVASCRIPT .... 112 6 1 MÉTODOS NUMÉRICOS NA SOLUÇÃO DE PROBLEMAS Segundo Quadros e Bortoli (2009), a maioria dos problemas da matemática é originária da necessidade de resolver situações da natureza. Numa primeira etapa, tem-se que obter um modelo matemático que representa de maneira conveniente um problema a ser analisado. Obtido o modelo matemático, procura-se encontrar a sua solução. Os mesmos autores citam que, quando se quer resolver um problema em engenharia, deve-se ter em mente o modelo que representa a situação física. Tal modelo é transformado em equações matemáticas (modelo matemático) que serão resolvidas por métodos analíticos ou por métodos numéricos. Como muitas situações não têm solução analítica ou a solução analítica é de difícil implementação computacional, os métodos numéricos se tornam a alternativa mais econômica em comparação à possibilidade de experimentação em laboratório. Este é um dos escopos do CÁLCULO NUMÉRICO. Outro objetivo é a programação computacional. Um exemplo importante é o cálculo de uma integral, frequentemente necessário em problemas da Engenharia Civil. Abaixo seguem duas integrais que têm solução analítica, ou seja, resolve-se substituindo a variável na solução por um número: Porém, existem integrais de funções que não têm solução analítica. Por exemplo, seja a integralde uma função f(x) no domínio a,b, que se quer resolver numericamente: Um dos métodos numéricos possíveis consiste em resolver a função f(x) em alguns pontos e somar as áreas sobre a curva da função, já que a integral nada mais é do que esse somatório: Figura: Área sobre a curva f(x) aproximada por um trapézio. Fonte: Quadros e Bortoli (2009). Do exemplo da figura anterior, pode-se entender também que todos os métodos numéricos são aproximações da solução analítica, com exatidão suficiente para seu propósito. Se dividirmos o intervalo a,b em vários intervalos menores, pode-se aumentar significativamente a exatidão do método. A integração numérica será estudada com mais detalhes posteriormente. 7 1.1 CÁLCULOS COMPUTACIONAIS: PROGRAMAÇÃO Obviamente, hoje em dia, a maioria dos cálculos numéricos de Engenharia não é feita com lápis e papel, já que lidamos com grande quantidade de cálculos e podemos utilizar os computadores. As calculadoras são usadas apenas para verificações simplificadas. Ocorre que existem muitos softwares com soluções prontas para vários tipos de situações e problemas, mas existem inúmeros problemas específicos que exigem cálculos também específicos. Podem-se citar aqui as variações regionais de materiais, mão-de-obra, climáticas, financeiras, etc. que exigem o desenvolvimento de soluções bem específicas. Citam-se ainda as pesquisas envolvendo novos materiais e novos processos de cálculo. Sendo assim, procura-se, nesta disciplina, fornecer as bases MATEMÁTICAS para o desenvolvimento de aplicativos computacionais. Sabe-se que as planilhas de cálculo (ex: Excel) possibilitam a resolução automatizada de inúmeras situações, por meio de fórmulas criadas para cada finalidade. Porém, casos que exigem maiores quantidades de cálculo e, especialmente, cálculos iterativos (“só param quando se atinge certa exatidão”) e cálculos em loop (“vai e volta”) são difíceis de realizar nas planilhas. Aplicativos de computador com entradas de dados mais “amigáveis” e possibilidades infinitas de cálculos são então a melhor solução na resolução dos problemas específicos de engenharia. Existem vários softwares (gratuitos e pagos) para a criação dos softwares, usando várias linguagens. Um bom software gratuito para este fim é o Lazarus (linguagem Pascal), que é livre e pode ser baixado do site: http://sourceforge.net/projects/lazarus/?source=dlp (acesso em 09/09/2013). O Lazarus permite a criação de aplicativos que rodam off-line (sem necessidade de internet), em qualquer computador. O aluno pode começar a estudar este software pelo conteúdo do Apêndice A desta apostila, que apresenta o passo-a-passo para a criação de um aplicativo simples (exercício 1.1 a seguir). Outros softwares/linguagens são FORTRAN e C++. Uma linguagem interessante para se aprender é a linguagem de programação dentro do Auto CAD, chamada AutoLisp (comando “vlisp”), sendo possível criar programas que geram e lêem desenhos dentro e arquivos dwg. Seja qual for a linguagem, o primeiro passo do aluno é aprender a lógica de programação, a montagem e resolução de algoritmos (sequências de cálculo), algo que independe da linguagem. O mesmo algoritmo pode ser implementado em qualquer linguagem, apenas se verificando como determinado comando é feito naquela linguagem. Atualmente, no mundo globalizado, é coerente criar os programas para execução num site da internet, on-line. Uma das possibilidades é usar a linguagem HTML combinada com o Javascript. Nesta disciplina, o aluno fica convidado a criar aplicações envolvendo as teorias estudadas em HTML/Javascript, sob a orientação do professor. Para tanto, o primeiro passo é aprender os conceitos básicos sobre a linguagem, sendo que o aluno pode fazer isso estudando, desde já, o Apêndice B desta apostila (exercício 1.2 a seguir). Aprendida uma linguagem, é relativamente fácil aprender qualquer outra, dependendo da futura necessidade profissional. O objetivo deste item é o primeiro contato do aluno com a programação. ************************************************************************************************ 8 TRABALHOS EM GRUPO: prezado aluno, todos os exercícios e criação de softwares desta disciplina serão em grupo. Cada grupo deve entregar um trabalho. Aproveitando a oportunidade, entenda que trabalhar em grupo significa se reunir com todos os integrantes e discutir cada exercício (ou parte do software) com todos, aproveitando as ideias dos colegas para criar a resposta final. Se fosse para cada integrante resolver um exercício (ou parte do software), o professor já faria a divisão desta maneira. Porém, se assim fosse, cada aluno só aprenderia a parte que fez. Aliás, nas avaliações, o conteúdo cobrado será de todas as partes de todos os trabalhos, e não partes deles. Aproveite que o curso é integral e que você tem a oportunidade de se reunir com o grupo, já que, profissionalmente, você não tem outro compromisso nesta fase da vida além de estudar (são raras as exceções de alunos que trabalham). Bom trabalho! EXERCÍCIOS 1.1) Criar o programa principal descrito no Apêndice A (área de triângulo). 1.2) Criar o programa principal descrito no Apêndice B (área de triângulo). ************************************************************************************************ 9 2 INTRODUÇÃO À TEORIA DE ERRO E ESTABILIDADE 2.1 FONTES DE ERROS A solução numérica dos problemas pode ser “contaminada” por alguns tipos de erro. ERRO DE TRUNCAMENTO Os erros de TRUNCAMENTO surgem exatamente quando se precisa parar a resolução quando o processo analítico seria infinito. Por exemplo, conforme Quadros e Bortoli (2009), podemos calcular numericamente o cos(x) pela resolução do polinômio abaixo (POLINÔMIO DE MCLAURIN), que vai até o infinito. Então é necessário truncar, parar num dado momento, gerando um erro. Obs: x tem que ser o ângulo em radianos. No caso acima, como conhecemos a solução exata, o erro “E” absoluto (neste caso, um erro de truncamento) é dado por: E = cos(x) – P(x) Em que: P(x) é uma função que aproxima a função exata. A tarefa é sempre, nas soluções, impor um erro máximo, uma tolerância “e”, e verificar se isso se satisfaz (E ≤ e), antes de continuar os cálculos. Obviamente, se estamos programando uma resolução numérica, é porque não existe solução analítica ou ela é inviável. Porém, sempre há uma maneira de estimar se a solução numérica converge para a solução exata. Isso é estudado de acordo com cada método numérico de solução, especificamente. Pode-se ainda calcular o erro relativo: )cos( )()cos( x xPx Er *100 (em %) É muito mais conclusivo se trabalhar com o erro relativo. Note que um erro de 1 km na medida da distância entre a Terra e a Lua é menos significativo que um erro de 1 cm na realização de uma cirurgia cerebral via robô. ************************************************************************************************ EXERCÍCIO 2.1) Na página do teu grupo (criada no exercício 1.2), crie um aplicativo que resolva o cosseno de um número segundo o explicado acima (polinômio de Mclaurin), digitando o código a seguir. Estude o código. (...) !6!4!2 1 642 xxx 10 ARQUIVO cosseno.html ARQUIVO funcoes-cos.js - Rode o programa off-line e verifique os resultados conforme abaixo: 11 - Faça testes aumentando cada vez mais o número de termos e perceba que, após muitos termos, o erro de truncamento não se altera. Ou seja, não seria necessário tanto esforço computacional para resolver o polinômiocom ótima precisão. Assim, você, programador, poderia limitar o número de termos em função de um erro aceitável, para reduzir esforço computacional. Para estudar mais aprofundadamente o cálculo do cosseno via polinômio de McLaurin, acesse: http://vtp.ifsp.edu.br/nev/Cosseno/cosseno.php? ************************************************************************************************* ERRO DE ARREDONDAMENTO Os erros de ARREDONDAMENTO surgem pelo fato de que as operações aritméticas quase nunca podem ser efetuadas com precisão completa, devido às aproximações dos números. Os tipos de arredondamento mais utilizados são, segundo Quadros e Bortoli (2009): - CORTE: as casas em excesso são simplesmente abandonadas; - MAIS PRÓXIMO: considerando-se “d” algarismos significativos, se o algarismo “d+1” for maior ou igual a 5, soma-se uma unidade ao algarismo de ordem “d”; caso contrário, o algarismo de ordem “d” permanece inalterado. Por exemplo, a solução exata da soma: Se uma máquina/aplicativo opera com 3 dígitos de precisão, considerando arredondamento por corte: Neste caso, o erro de arredondamento por corte é igual a 2-1,999=0,001. Na computação atual, as máquinas trabalham com muitos dígitos de precisão e o erro de arredondamento é automaticamente minimizado. A questão da precisão fica somente para os resultados numéricos mostrados pelo programa. Por exemplo, no programa anterior, veja que os valores mostrados do cosseno exato e do polinômio de McLaurin estão mostrados com todas as casas decimais padrão. O programador pode alterar essa precisão mostrada, de acordo com a aplicação e público-alvo. ************************************************************************************************* 12 EXERCÍCIO 2.2) Altere o aplicativo do exercício 2.1 e defina a precisão dos resultados para: Cosseno analítico precisão de mais próximo, 5 casas decimais Cosseno numérico precisão de mais próximo, 5 casas decimais Erro precisão de corte, 6 casas decimais As funções que podem ser usadas são: n.toPrecision(5) // “n” com precisão de mais próximo em 5 algarismos (total de algarismos, incluindo decimais) n.toFixed(5) // “n” com precisão de mais próximo em 5 casas decimais Math.floor(n*Math.pow(10,5))/Math.pow(10,5); //“n” com precisão de corte em 5 casas decimais ************************************************************************************************* 2.2 PRECISÃO E EXATIDÃO Torna-se aqui necessário distinguir precisão de exatidão. Precisão é relacionada à quantidade de algarismos significativos. Exatidão é a medida da perfeição do resultado. A exatidão depende da precisão da máquina/aplicativo e do método utilizado para a obtenção do resultado. Veja o exemplo de Quadros e Bortoli (2009), envolvendo o número que vale 3,1415926: 2.3 ESTABILIDADE Um método numérico diz-se instável se a acumulação de pequenos erros durante o cálculo pode ter grande influência na precisão dos resultados. Um método estável produz bons resultados mesmo com o acúmulo de pequenos erros. Chamam-se os algoritmos (sequência de solução) instáveis também de algoritmos mal condicionados. Cherri et. al. (????) exemplifica o seguinte problema instável: supondo-se que as operações abaixo sejam processadas em uma máquina ou aplicativo com 4 dígitos significativos: 13 ALGORITMO A) ALGORITMO B) Os dois resultados são diferentes, quando não deveriam ser, pois a adição é uma operação distributiva e, obviamente, o resultado é x2. A causa dessa diferença foi um arredondamento feito na adição (x2 + x1), cujo resultado tem 8 dígitos. Como a máquina, numa hipótese deste exemplo, só armazena 4 dígitos, os menos significativos foram desprezados. Por isso é importante somente arredondar os números no momento em que for apresentado o resultado final, deixando o programa trabalhar nos cálculos internos com todas as casas decimais possíveis. Concluindo, ao se utilizar softwares de programação ou máquinas calculadoras, deve-se ficar atento aos erros que podem acontecer, sejam por truncamento ou por arredondamento, que podem levar a um algoritmo instável. 14 3 SISTEMAS DE EQUAÇÕES LINEARES 3.1 SISTEMAS DE EQUAÇÕES LINEARES NA ENGENHARIA CIVIL Podemos aplicar o cálculo numérico, por exemplo, à análise de estruturas de Engenharia Civil, no cálculo dos esforços. Um dos métodos numéricos existentes é o Método dos Elementos Finitos (MEF)1, que reduz o cálculo de esforços à resolução de um sistema de equações lineares, com a matriz (chamada matriz de rigidez) vezes o vetor de deslocamentos igual ao vetor de forças ( GG PK * ). A estrutura abaixo, por exemplo, quando calculada via MEF, resume-se a resolver o sistema linear 4x4 exposto na sequência (para rigidez E*I = 1 kN*m2). 0 0,4 0 4,0 '3 3 '2 2 * 8026 024612 2680 612024 v v v v (FORMA MATRICIAL) Em forma de equações, o sistema fica: 24*v2 -12*v3 + 6*v3‟ = 0,4 (1) 8*v2‟ -6*v3 + 2*v3‟ = 0 (2) -12*v2 -6*v2‟ +24*v3 = 0,4 (3) 6*v2 +2*v2‟ +8*v3‟ = 0 (4) Resolvendo o sistema, os deslocamentos, ficam: v2 = -0,0667 m (deslocamento vertical do nó 2) v2‟= -0,0667 rad (giro do nó 2) v3 = -0,0667 m (deslocamento vertical do nó 3) v3‟ = 0,0667 rad (giro do nó 3) Com os parâmetros calculados, é possível, posteriormente, calcular os esforços da estrutura. O que o aluno tem que entender aqui não é o MEF em si, mas que o MEF utiliza a resolução de um sistema linear de equações, o que se vai estudar a seguir. Cumpre destacar que há muitos outros métodos em Engenharia Civil, em suas diversas subáreas, que precisam da resolução de sistemas lineares. 1 O MEF será estudado em outra(s) disciplina(s) ou na pós-graduação. 15 3.2 SOLUÇÃO NUMÉRICA COMPUTACIONAL DE SISTEMAS LINEARES DE EQUAÇÕES Um sistema linear da forma AX = B pode ser resolvido por meio de métodos diretos ou de métodos iterativos. Um método é direto quando a solução exata é obtida realizando-se um número finito de operações aritméticas. Um método é iterativo quando a solução é obtida como limite de uma sequência de aproximações sucessivas, ou seja, deve-se truncar, parar a solução após certo número de iterações. 3.2.1 MÉTODOS DIRETOS PARA A SOLUÇÃO DE SISTEMAS LINEARES Considere um sistema triangular da forma: Ele pode ser resolvido por retro substituição, ou seja, resolve-se a última linha, substituindo o resultado na penúltima, e assim por diante. Como regra geral, podem-se assim obter as soluções: ; A seguir apresentam-se os algoritmos para resolução de sistemas triangulares superiores (exemplo acima, onde se começa a resolução pela última linha) e inferiores (contrário, onde a primeira linha só tem 1 termo e a resolução começa por ela): RESOLUÇÃO DE SISTEMAS TRIANGULARES SUPERIORES 1: x(n) = b(n)/a(n,n) 2: Para i=n-1 até i=1 faça //a partir da penúltima linha, até a primeira 3: soma=0 4: Para j=i+1 até j=n faça 5: soma=soma + a(i,j)*x(j) 6: Fim do laço 7: x(i) = (b(i)-soma) / a(i,i) 8: Fim do laço 16 RESOLUÇÃO DE SISTEMAS TRIANGULARES INFERIORES 1: x(1) = b(1)/a(1,1) 2: Para i=2 até i=n faça //a partir da segunda linha, até a última 3:soma=0 4: Para j=1 até j=i-1 faça 5: soma=soma + a(i,j)*x(j) 6: Fim do laço 7: x(i) = (b(i)-soma)/a(i,i) 8: Fim do laço Daqui para frente o aluno irá se deparar com algoritmos, que são as sequências de cálculo de uma determinada tarefa. Quando se quer entender o que um algoritmo faz, pode-se calculá-lo passo-a-passo. Segue um exemplo usando o algoritmo para resolução de SISTEMAS TRIANGULARES SUPERIORES apresentado: n=4 x4=b4/a44=2/2=1 i=n-1=3 soma=0 j=i+1=4 (=n, portanto último e único passo) soma=soma+a34*x4=0+(-5)*1=-5 x3=(b3-soma)/a33=(3-(-5))/4=8/4=2 i=i-1=2 soma=0 j=i+1=3 soma=soma+a23*x3=0+1*2=2 j=j+1=4 (=n, portanto último passo) soma=soma+a24*x4=2+(-2)*1=0 x2=(b2-soma)/a22=((-1)-0)/1=-1 i=i-1=1 (=n, portanto último passo) soma=0 j=i+1=2 soma=soma+a12*x2=0+4*(-1)=-4 j=j+1=3 soma=soma+a13*x3=(-4)+(-5)*2=-14 j=j+1=4 (=n, portanto último passo) soma=soma+a14*x4=(-14)+1*1=-13 x1=(b1-soma)/a11=((-10)-(-13))/3=3/3=1 Esta resolução passo-a-passo verificando a variação do valor de cada variável ao longo da execução do algoritmo é chamada DEPURAÇÃO (“DEBUG”, em inglês, referindo-se ao ato de procurar erros, “bugs”). 1 2 1- 1 x 17 Normalmente os softwares para criar softwares tem uma função debug, em que apresenta uma tabela mostrando a variação dos valores das variáveis numa execução passo-a-passo. “PSEUDOCÓDIGO” Se o algoritmo for apresentado genericamente, como o exemplo anterior, ou seja, com comandos “para”, “fim”, etc., em português, e não nos comandos exatos de uma linguagem de programação específica, é chamado PSEDOCÓDIGO. O que ocorre é que cada linguagem usa determinado nome para cada comando, sendo que o mesmo algoritmo tem que ser escrito com os comandos da linguagem de programação que será implementada. Por exemplo, uma iteração do tipo “para i=1 até n” é criada como “do i=1 to n” em linguagem FORTRAN e “for (i=1; i<n; i++)” em Javascript. Uma forma equivalente e, algumas vezes, complementar ao pseudocódigo, utilizada para se representar um algoritmo é o diagrama de fluxos (fluxograma). O Fluxograma utiliza figuras geométricas para representar um algoritmo. Algumas das formas básicas são: Início ou fim do algortimo condicional procedimentos, contas, etc. Leitura de dados saída de dados (impressão ou em tela) Segue um exemplo simples de algoritmo criado com fluxograma: 18 ************************************************************************************************* EXERCÍCIOS 3.1) Resolver, em papel, usando o algoritmo apresentado, o seguinte sistema linear triangular superior de ordem 4 (apresentar o passo-a-passo com todas as variáveis; pode ser na forma de tabela de depuração, mas não é obrigatório): 3x1+2x2-4x3+3x4=-5 6x2-x3+x4=14 2x3+x4=11 5x4=10 3.2) Resolver, em papel, usando o algoritmo apresentado, o seguinte sistema linear triangular inferior de ordem 4 (apresentar o passo-a-passo com todas as variáveis): 2x1=8 5x1+2x2=16 x1-x2+x3=5 3x1+5x2-2x3+2x4=19 3.3) Crie um programa em HTML/Lazarus, na página do GRUPO, para resolver um sistema linear triangular qualquer superior de ordem n (n x n) qualquer. Abaixo segue um algoritmo possível para a página HTML e arquivo .js com algoritmos para a leitura e impressão dos dados (sem os cálculos). O aluno deve implementar o algoritmo de resolução descrito. ARQUIVO sistemasuperior.html <html> <meta http-equiv="Content-Type" content="text/html;charset=utf-8" /> <title>Sistema Linear Triangular Superior</title> <script src='funcoes-sistemasuperior.js'></script> SISTEMA LINEAR TRIANGULAR SUPERIOR <br><br> Ordem do sistema (n): <input type="number" id="ordem" value="4"> <br><br> <button onclick="gerar()">GERAR</button> <br><br> Obs: Matriz COMPLETA (a última coluna - cinza - contém as constantes, os valores após a igualdade). <br><br> <ins id="tabela"></ins> <!-- TABELA IMPRESSA, DE ACORDO COM A FUNÇÃO gerar() no arquivo .js --> <br> <button onclick="calcular()">CALCULAR</button> <br><br> <ins id="resultado"></ins> <!-- IMPRESSÃO DOS RESULTADOS, DE ACORDO COM A FUNÇÃO calcular() no arquivo .js --> <br><br> </html> 19 ARQUIVO funções-sistemasuperior.js function gerar(){ n=document.getElementById("ordem").value; //lê n n=Number(n); //para definir que n é número, não texto col=n+1; //número de colunas, porque a última coluna são os valores após a igualdade //abaixo, html dentro do javascript, para criar a tabela dinamicamente (nxn) html=""; html+='<table border>'; for (i=1;i<=n;i++){ html+='<tr>'; for(j=1;j<=n;j++){ ntermo=i+','+j; conteudo='<input type="number" id='+ntermo+' value="0"/>'; html+='<td>'+conteudo+'</td>'; } ntermo=i+','+col; conteudo='<input type="number" style="background-color: #C0C0C0" id='+ntermo+' value="0"/>'; html+='<td>'+conteudo+'</td>'; html+='</tr>'; } html+='</table>'; document.getElementById('tabela').innerHTML=html; } function calcular(){ n=document.getElementById("ordem").value; n=Number(n); //ler os termos aij a=[ ]; //criar o vetor for (i=1;i<=n;i++){ for(j=1;j<=n;j++){ ntermo=i+','+j; a[ntermo]=document.getElementById(""+ntermo+"").value; a[ntermo]=Number(a[ntermo]); } } //ler os termos bi col=n+1; b=[ ]; //criar o vetor for (i=1;i<=n;i++){ ntermo=i+','+col; b[i]=document.getElementById(""+ntermo+"").value; b[i]=Number(b[i]); } sistemasuperior(n,b,a); //chama a função que calcula o sistema triangular superior //imprimir os resultados xi res=''; for (i=1;i<=n;i++){ res=res+'\r\nx'+i+'\r\n='+x[i]+'<br>' } document.getElementById('resultado').innerHTML = res; } function sistemasuperior(n,b,a){ //AQUI ENTRA O ALGORITMO DOS ALUNOS //entra n, b e a. Sai x. IMPLEMENTE O ALGORITMO AQUI. (...) } 20 O exemplo cujo passo-a-passo do algoritmo foi mostrado anteriormente está a seguir, no formato que deve ficar o programa: 3.4 (EXTRA)) Duplique o programa anterior e o altere para resolver um sistema linear triangular qualquer inferior de ordem n. Um exemplo adaptado a partir do anterior, invertendo as equações, fica: 2x1=2 -5x1+4x2=3 -2x1+x2+x3=-1 x1-5x2+4x3+3x4=-10 ************************************************************************************************* 1 1- 2 1 x 21 As soluções dos sistemas lineares inferior e superior são necessárias no MÉTODO DA ELIMINAÇÃO DE GAUSS. MÉTODO DA ELIMINAÇÃO DE GAUSS (OU MÉTODO DO ESCALONAMENTO) Este método consiste em construir, a partir de um sistema de equações lineares qualquer, um sistema triangular A′X = B′ equivalente que pode ser resolvido por retro substituição, com os algoritmos estudados anteriormente. Quadros e Bortoli (2009) citam que dois sistemas lineares de dimensão n x n são equivalentes desde que os seus conjuntos de soluções sejam os mesmos. Teoremas de álgebra linear mostram que, quando certas operações são aplicadas a um dado sistema, a solução não muda. Citam os autores que as seguintes operações, quando aplicadas a um sistema linear, resultam num sistema equivalente: 1. Mudança de ordem de duas equações. 2. Multiplicação de uma equação por uma constante não nula. 3. Adição de um múltiplo de uma equação a outra equação. Quadros e Bortoli (2009) descrevemo método de eliminação de Gauss para um sistema de ordem 3, sendo que o mesmo processo pode ser aplicado a sistemas de qualquer ordem. Apresentam o seguinte sistema: Considere-se somente a matriz dos coeficientes aij aumentada contendo os termos independentes bi: O objetivo é obter um sistema triangular da forma abaixo, sendo que os elementos da segunda até a última linha (no caso, a terceira), obviamente, não serão os mesmos da matriz original (por isso estão indicados com o índice linha): O processo consiste em fazer uma combinação linear para zerar os elementos da primeira coluna abaixo do pivô (primeiro elemento – da esquerda para a direita - não-nulo de cada linha; por isso se diz estar fazendo um PIVOTEAMENTO). Para 22 tanto, segundo Cantão (????), para cada linha “i”, a partir da segunda, subtrai-se a primeira linha multiplicada pelo fator: Quadros e Bortoli (2009) apresentam o seguinte exemplo: A matriz aumentada para este sistema é (o pivô está sublinhado): m21 = a21/a11 = 2/1 = 2 m31 = a31/a11 = 4/1 = 4 m41 = a41/a11 = -3/1 = -3 a‟21 = a21 – m21*a11 = 2-2*1 = 0 (ok) a‟22 = a22 – m21*a12 = 0-2*2 = -4 a‟23 = a23 – m21*a13 = 4-2*1 = 2 a‟24 = a24 – m21*a14 = 3-2*4 = -5 b‟2 = b2 – m21*b1 = 28-2*13 = 2 a‟31 = 0 a‟32 = a32 – m31*a12 = 2-4*2 = -6 a‟33 = a33 – m31*a13 = 2-4*1 = -2 a‟34 = a34 – m31*a14 = 1-4*4 = -15 b‟3 = b3 – m31*b1 = 20-4*13 = -32 a‟41 = 0 a‟42 = a42 – m41*a12 = 1-(-3)*2 = 7 a‟43 = a43 – m41*a13 = 3-(-3)*1 = 6 a‟44 = a44 – m41*a14 = 2-(-3)*4 = 14 b‟4 = b4 – m41*b1 = 6-(-3)*13 = 45 Portanto, o resultado, após a eliminação, para a 2ª matriz, é (já informando os parâmetros pivô e “m” para a próxima etapa): 23 A terceira matriz fica: m32 = a32/a22 = -6/-4 = 1,5 m42 = a42/a22 = 7/-4 = -1,75 a‟31 = 0 a‟32 = 0 a‟33 = a33 – m32*a23 = (-2)-1,5*2 = -5 a‟34 = a34 – m32*a24 = (-15)-1,5*(-5) = -7,5 b‟3 = b3 – m32*b2 = (-32)-1,5*2 = -35 a‟41 = 0 a‟42 = 0 a‟43 = a43 – m42*a23 = 6-(-1,75)*2 = 9,5 a‟44 = a44 – m42*a24 = 14-(-1,75)*(-5) = 5,25 b‟4 = b4 – m42*b2 = 45-(-1,75)*2 = 48,5 A quarta e última matriz fica: m43 = a43/a33 = 9,5/-5 = -1,9 a‟41 = 0 a‟42 = 0 a‟43 = 0 a‟44 = a44 – m43*a34 = 5,25-(-1,9)*(-7,5) = -9 b‟4 = b4 – m43*b3 = 48,5-(-1,9)*(-35) = -18 Observando que todos os elementos da diagonal são não nulos, a solução, calculada conforme exposto anteriormente (sistema triangular superior), fica: Para n = 4; i = 4: 2 9 18 44 4 4 a b x Para n = 4; i = 3: 24 4 5 )2*5,7(35 * 33 4 4 33 3 a xab x j jj Para n = 4; i = 2: 1 4 ))2*)5((4*2(2 * 22 4 3 22 2 a xab x j jj Para n = 4; i = 1: 3 1 ))2*4()4*1())1(*2((13 * 11 4 2 11 1 a xab x j jj A montagem do algoritmo, seguindo os conceitos expostos, vem na sequência. Quadros e Bortoli (2009) descrevem o método de eliminação de Gauss para um sistema de ordem 3, sendo que o mesmo processo pode ser aplicado a sistemas de qualquer ordem. Partindo-se da matriz aumentada, na primeira combinação, o objetivo é: Como visto no exemplo, a combinação linear (considerando os critérios 1 a 3 dispostos no começo do capítulo) que chega ao resultado da matriz anterior (zerando a coluna 1) é, segundo Cantão (????), para cada linha “i”, a partir da segunda, subtrair a primeira linha multiplicada pelo fator: Cada elemento da linha i modificada fica (o expoente “(2)” é para indicar que os elementos são de uma 2ª matriz): Para i=2,n Para j=1,n (sendo que, para j=1, aij (2) dará zero) Sendo assim, os novos valores para o exemplo de 3 linhas são dados por: 25 linha 2 linha 3 a'21 = a21 – m21a11 = 0 a'31 = a31 – m31a11 = 0 Repete-se o processo para zerar os elementos inferiores da segunda coluna. Então, zerando os elementos da coluna k, para cada linha (i > k), que tem o pivô, o multiplicador “m” o cálculo de dos elementos assumem a forma genérica: natéki 1 //para cada linha depois da anterior, que tem o pivô KK IK IK a a m (atenção: os coeficientes “a” são os modificados da etapa anterior) KIK K I K I bmbb 1 natékj //para cada coluna, a partir da k+1 KJIK K IJ K IJ amaa 1 Para zerar os elementos da segunda coluna (k=2), os elementos da matriz ficam: linha 3 a'‟32 = a‟32 – m32a‟22 = 0 O processo é conhecido também como ESCALONAMENTO. Quadros e Bortoli (2009) lembram que o sistema final não pode ter elementos nulos na diagonal, ou seja, os elementos a11, a22, a33 etc. não podem ser nulos, o que resultaria num sistema sem solução, a princípio, por este método. Conforme Quadros e Bortoli (2009), o método de Gauss irá falhar quando um pivô for nulo, pois, neste caso, não seria possível calcular os multiplicadores “m” Este sério problema pode ser evitado pelo uso da estratégia de pivoteamento PARCIAL. Esta consiste em trocar linhas (ou colunas) de forma a ter sempre o pivô não-nulo. A técnica de pivoteamento parcial tem o seguinte esquema de escolha do pivô: 1º Pivô = elemento de maior valor absoluto (módulo) na coluna 1; 2º Pivô = elemento de maior valor absoluto na coluna 2; E assim sucessivamente. A seguir apresenta-se um exemplo com pivoteamento parcial: 26 A matriz aumentada para este sistema é: Fazendo-se o pivoteamento parcial, ou seja, colocando a linha com o maior pivô como a primeira: 4584 3742 3309 Fazendo uma operação para zerar os termos abaixo da linha do pivô, conforme explicado anteriormente: m21 = a21/a11 = 2/9 = 0,222 m31 = a31/a11 = 4/9 = 0,444 a‟21 = 0 a‟22 = a22 – m21*a12 = (-4)-0,222*0 = -4 a‟23 = a23 – m21*a13 = 7-0,222*(-3) = 7,666 b‟2 = b2 – m21*b1 = 3-0,222*3 = 2,334 a‟31 = 0 a‟32 = a32 – m31*a12 = (-8)-0,444*0 = -8 a‟33 = a33 – m31*a13 = 5-0,444*(-3) = 6,332 b‟3 = b3 – m31*b1 = (-4)-0,444*3 = -5,332 332,5332,680 334,2666,740 3309 Fazendo-se novamente o pivoteamento parcial, ou seja, colocando a linha com o maior pivô (em módulo) como a próxima linha de resolução: 334,2666,740 332,5332,680 3309 Para o segundo pivô: m32 = a32/a22 = -4/-8 = 0,5 27 a‟31 = 0 a‟32 = 0 a‟33 = a33 – m32*a23 = 7,666-0,5*6,332 = 4,5 b‟3 = b3 – m32*b2 = 2,334-0,5*(-5,332)= 5 55,400 332,5332,680 3309 E, finalmente, resolvendo o sistema triangular superior: x1 = 0,7037; x2 = 1,5462; x3 = 1,1111 A seguir apresenta-se um algoritmo para a resolução de sistema de equações lineares pelo Método de Eliminação de Gauss com pivoteamento parcial: Dado n, Anxn e bnx1 1: Para k=1 até k=n-1 faça //para cada coluna (a última de “a” e a coluna “b” não entram) //Encontrar o pivô (o maior em módulo) 2: w=|a(k,k)| //w é o pivô 3: r=k //r é o nº da linha que está o pivô 4: Para j=k até j=n faça //para cada linha menos as anterioresa “k”, já escalonadas 5: Se |a(j,k)|>w então 6: w=|a(j,k)| 7: r=j 8: Fim do condicional 9: Fim do laço 10: Se w=0 então 11: „Todos os pivôs são nulos.‟ PARE. 12: Senão //Fim encontrar o pivô (o maior em módulo) //Trocar a linha k com a linha r 13: Para c=1 até c=n faça //coluna 1 a n 14: atemp=a(k,c) //atemp é uma variável temporária 15: a(k,c) = a(r,c) 16: a(r,c) = atemp 17: Fim do laço 18: btemp = b(k) //btemp é uma variável temporária 19: b(k) = b(r) 20: b(r) = btemp 21: Fim do condicional //Fim Trocar a linha k com a linha r //Fazer o escalonamento 22: Para i=k+1 até i=n faça //para cada linha depois da que tem o pivô 23: m(i,k)=a(i,k)/a(k,k) 24: b(i)= b(i) - m(i,k)* b(k) 25: a(i,k)=0 //vai dar zero, mas é bom colocar para evitar erro de precisão 26: Para j=k+1 até j=n faça //para cada coluna depois da que tem o pivô 27: a(i,j)= a(i,j) - m(i,k)* a(k,j) 28: Fim do laço 29: Fim do laço 30: Fim do laço //Fim fazer o escalonamento (“a” já está escalonado, ou seja, já é triangular superior) 31: Execute o algoritmo de solução de sistemas triangulares superiores com a,b,n. 28 Para estudar mais aprofundadamente o Método da Eliminação de Gauss com pivoteamento parcial, acesse: http://vtp.ifsp.edu.br/nev/Sistema-gauss/sistemagauss.php? ************************************************************************************************* 29 EXERCÍCIOS 3.5) Resolver, em papel, o seguinte sistema linear pelo Método de Gauss, com pivoteamento parcial: 3.6) Resolver, em papel, o exercício anterior usando o algoritmo apresentado, apresentando o passo a passo. 3.7) Duplique o programa do exercício 3.3 (sistema triangular superior) e modifique-o para resolver sistemas pelo algoritmo de eliminação de Gauss com pivoteamento parcial, que acabou de ser apresentado. Para implementar o algoritmo, seguem algumas informações sobre a linguagem Javascript, na ordem do algoritmo proposto: Módulo de um número n: Math.abs(n) Condicional: if (a<0) { //qualquer código } Obs: use os operadores lógicos a seguir: == igual a //ATENÇÃO, NA COMPARAÇÃO COLOQUE 2 IGUAIS != diferente >, <, >=, <= maior, menor, maior ou igual e menor ou igual PARE (parar o programa): return; Condicional com “SENÃO”: if (a==0) { //código } else { //código } ************************************************************************************************* 30 MÉTODO DA FATORAÇÃO LU Lembre-se que, por exemplo, o Método dos Elementos Finitos utiliza a resolução de um sistema linear para encontrar os deslocamentos nodais () de uma estrutura. De acordo com o MEF, o sistema é formado pela matriz de rigidez (K) da estrutura e pelos vetores de cargas nodais equivalentes (b), formando o sistema linear K*=b. Acontece que as estruturas civis podem ser tão grandes que a ordem do sistema chega à casa dos milhares. Além disso, para a mesma estrutura (mesma matriz K), devem ser consideradas dezenas de combinações de ações (vetor b). Sendo assim, pelo escalonamento de Gauss estudado, devem-se resolver dezenas de sistemas lineares da ordem de milhares, escalonando-o a cada resolução. Computacionalmente, isto se torna bastante demorado. Então, será a seguir estudado um método que permite escalonar apenas a matriz de rigidez da estrutura (K) (uma única vez) e, para cada vetor de cargas nodais equivalentes (b), resolver 2 sistemas lineares triangulares (sem qualquer escalonamento). Sendo assim, o método da fatoração LU é um método derivado do Método da Eliminação de Gauss, que apresenta a vantagem de que o escalonamento é feito independentemente do vetor dos termos independentes bi. Então se pode escalonar a matriz dos coeficientes e salvá-la na memória, resolvendo um sistema equivalente para qualquer vetor de termos independentes bi. No processo de Gauss, toda vez que se muda o vetor bi, deve-se fazer novamente todo o escalonamento. De acordo com Ruggiero (1996), neste método, para resolver o sistema linear Ax=b, se pudermos realizar a fatoração A=CD, o sistema fica (CD)x=b. Se y=Dx, então resolver o sistema é equivalente a resolver Cy=b e, depois, Dx=y. Cantão (????) cita que LU é um termo que vem do inglês lower e upper, já que o método se baseia na decomposição da matriz A na forma A = L*U, onde L é uma matriz triangular inferior (todos os elementos acima da diagonal são nulos) com elementos da diagonal principal iguais a 1 e U é uma matriz superior (todos os elementos abaixo da diagonal são nulos) qualquer. Primeiro vamos teorizar o método, usando um sistema linear literal, conforme Ruggiero (1996): Considerando o escalonamento somente da matriz dos coeficientes A, sem o vetor das constantes b, porque queremos fatorar somente A, conforme anteriormente explicado, os multiplicadores da etapa 1 são: E, normalmente, os coeficientes para gerar a matriz A(1) são calculados como (exatamente da mesma maneira que feito no Método de Gauss): 31 Percebe-se que os mesmos resultados são conseguidos, para a matriz A(1) da 1ª etapa, multiplicando A(0) por M(0) conforme abaixo: Para a 2ª etapa, supondo (exatamente da mesma maneira que feito no Método de Gauss): Percebe-se que os mesmos resultados são conseguidos, para a matriz A(2) da 2ª etapa, multiplicando A(1) por M(1) conforme abaixo: Em resumo: A(1) = M(0)A(0) A(2) = M(1)A(1) = M(1)M(0)A(0) Então: A(0) = (M(1))(-1)(M(0))(-1)A(2) = A (2) 32 Então, sendo A a matriz original (ou A(0)), encontra-se a fatoração LU procurada: Então, resolver o sistema Ax=b é equivalente a resolver Ly=b e, depois, Ux=y. Note que a dedução não considerou o pivoteamento parcial. Cantão (????) apresenta um exemplo: A*x = b Para a matriz A, o escalonamento é feito normalmente como apresentado no Método de Eliminação de Gauss com pivoteamento normal (não parcial): 1ª Etapa: m21 = a21/a11 = 1/3 m31 = a31/a11 = 4/3 a‟21 = 0 a‟22 = a22 – m21*a12 = 1-1/3*2 = 1/3 a‟23 = a23 – m21*a13 = 2-1/3*4 = 2/3 a‟31 = 0 a‟32 = a32 – m31*a12 = 3-4/3*2 = 1/3 a‟33 = a33 – m31*a13 = 2-4/3*4 = -10/3 33 2ª Etapa: m32 = a32/a22 = (1/3)/(1/3) = 1 a‟31 = 0 a‟32 = 0 a‟33 = a33 – m32*a23 = (-10/3)-1*(2/3) = -4 A escalonada = U = Montando a partir dos multiplicadores calculados L = = Resolvendo Ly=b: Na resolução anterior, lembre-se que o expoente T significa que o vetor y é transposto, ou seja, que, na verdade, os valores estão em coluna, não em linha. Resolvendo Ux=y: Exposto o método, obviamente é interessante resolver o Método considerando o pivoteamento parcial. Ruggiero (1996) diz que o impacto do pivoteamento parcial é somente fazer a mesma troca de linhas feita na matriz A, no escalonamento, no vetor b quando da resolução do sistema Ly=b. E esta troca de linhas também deve ser realizada na matriz L, mas somente em sua parte inferior, sendo que a diagonal se mantém unitária e a parte superior se mantém zerada. Na prática, o que se faz é, durante o escalonamento, criar um vetor P com a ordem das linhas e e, a cada troca, trocar a ordem deste vetor. Ao final do escalonamento, ter- se-á a troca que deve ser feita no vetor b. 34 RESUMO DO PROCESSO LU COM PIVOTEAMENTO PARCIAL (resolver o sistema linear Ax=b): - escalonar, via Gauss com pivoteamentoparcial, A até chegar a U. Portanto, U é a matriz A escalonada (triangular superior) - a cada etapa (“k”) do pivoteamento parcial, faz-se uma troca de linhas, armazenando as trocas num vetor P, que é alterado a cada etapa - portanto, ao final do escalonamento, além de U, obtém L (matriz triangular inferior com os coeficientes multiplicadores m, já com as trocas de linha na parte inferior) e P com a troca final de linhas - Ly=Pb obtém y (sistema linear triangular inferior) - Ux=y obtém x (que é o objetivo; sistema linear triangular superior) Ruggiero (1996) apresenta um exemplo: P = (1,2,3) (ordem inicial) Fazendo o pivoteamento (1ª etapa), o maior pivô está na linha 3 e, portanto, devemos trocá-la com a 1ª linha: e P=(3,2,1) (troca com relação a P anterior) Ao efetuar a eliminação de Gauss, Ruggiero (1996) usa as posições zeradas para guardar os multiplicadores m. Isso é comum em algoritmos para economizar memória computacional. Dessa maneira, efetuando a eliminação na matriz A‟(0): Fazendo o pivoteamento (2ª etapa), o maior pivô está na linha 3 e, portanto, devemos trocá-la com a 2ª linha: e P=(3,1,2) (troca com relação a P anterior) 35 Efetuando a eliminação na matriz A‟(1): E, portanto: P=(3,1,2) Pb é o mesmo que trocar as linhas de b segundo P. Sendo assim, conforme P, a primeira linha de Pb será a linha 3 de b, a segunda linha será a linha 1 e a terceira linha será a linha 2: bT = (9, 3, -2) (Pb)T = (-2,9,3) O algoritmo para a resolução por fatoração LU está descrito a seguir. Estão numeradas apenas as linhas acrescentadas ou que mudam em relação ao algoritmo apresentado para resolução via Gauss com pivoteamento parcial. 36 Dado n, Anxn e bnx1 //criação do vetor p inicial (com a ordem das linhas) – 1,2,3,...,n 1: Para i=1 até i=n faça 2: p(i)=i 3: Fim do laço //fim da criação do vetor p inicial (com a ordem das linhas) – 1,2,3,...,n //escalonamento da matriz A Para k=1 até k=n-1 faça //para cada coluna (a última de “a” não entra) //Encontrar o pivô (o maior em módulo) w=|a(k,k)| //w é o pivô r=k //r é o nº da linha que está o pivô Para j=k até j=n faça //para cada linha menos as anteriores a “k”, já escalonadas Se |a(j,k)|>w então w=|a(j,k)| r=j Fim do condicional Fim do laço Se w=0 então „Todos os pivôs são nulos.‟ PARE. Senão //Fim encontrar o pivô (o maior em módulo) //Trocar a linha k com a linha r Para c=1 até c=n faça //coluna 1 a n atemp=a(k,c) //atemp é uma variável temporária a(k,c) = a(r,c) a(r,c) = atemp Fim do laço //trocar as linhas do vetor p (não do b, como é no Gauss) 4: ptemp = p(k) //ptemp é uma variável temporária 5: p(k) = p(r) 6: p(r) = ptemp Fim do condicional //trocar as linhas do vetor p //Fim Trocar a linha k com a linha r //Fazer o escalonamento Para i=k+1 até i=n faça //para cada linha depois da que tem o pivô m(i,k)=a(i,k)/a(k,k) 7: //tirada esta linha com escalonamento de b 8: a(i,k)= m(i,k) //armazenando os multiplicadores na própria matriz A escalonada, no lugar de 0 Para j=k+1 até j=n faça //para cada coluna depois da que tem o pivô a(i,j)= a(i,j) - m(i,k)* a(k,j) Fim do laço Fim do laço Fim do laço //fim escalonamento da matriz A (sendo sua parte inferior, a matriz L e sua parte superior a matriz U 9: //aqui não mais executa diretamente a resolução de um sistema triangular superior //troca das linhas de b, conforme o vetor p que tem as trocas de linhas feitas em A. Estou chamando “Pb” de “c” 10: Para i=1 até i=n faça 11: r=p(i) 12: c(i)=b(r) 13: Fim do laço //fim troca das linhas de b, conforme o vetor p que tem as trocas de linhas feitas em A //resolver o sistema triangular inferior Ly=c – como L está em A, mas não tem a diagonal =1, só alteramos o algoritmo de sistema triangular inferior tirando a divisão por a(i,i) na última linha (linha 20 a seguir), porque divisão por 1 não é necessária. Alteram-se também as variáveis 14: y(1) = c(1) //não divide por a(1,1), que é 1 15: Para i=2 até i=n faça //a partir da segunda linha, até a última 16: soma=0 37 17: Para j=1 até j=i-1 faça 18: soma=soma + a(i,j)*y(j) 19: Fim do laço 20: y(i) = (c(i)-soma) //aqui tirou a divisão por a(i,i), que é 1 21: Fim do laço //fim resolver o sistema triangular inferior Ly=c //resolver Ux=y (sistema triangular superior) 22: Execute o algoritmo de solução de sistemas triangulares superiores com a,y,n. // veja que pode ser usada a mesma função que resolve sistema triangular superior, só que entra com a,y,n e não a,b,n. “A” já é “U” na parte superior e o algoritmo só considera mesmo a parte superior ************************************************************************************************* Para estudar mais aprofundadamente o Método de LU, acesse: http://vtp.ifsp.edu.br/nev/Sistema-lu/sistemalu.php? ************************************************************************************************* 38 EXERCÍCIOS: 3.8) Resolver, em papel, o seguinte sistema linear (o mesmo do exercício 3.5) pela fatoração LU com pivoteamento normal (não parcial): 3.9) Resolver, em papel, o sistema do exercício anterior pela fatoração LU com pivoteamento parcial. 3.10) Resolver, em papel, o exercício anterior usando o algoritmo apresentado, apresentando o passo a passo. 3.11) Usando como base os programas anteriormente criados, crie um software para resolver sistemas pelo algoritmo de fatoração LU com pivoteamento parcial, que acabou de ser apresentado. O software deve mostrar as matrizes completas L|c (sendo que L é a parte inferior de A, exceto a diagonal com elementos iguais a 1) e U|y (sendo que U é a parte superior de A, inclusive com a diagonal; a parte inferior é zerada) conforme o exemplo a seguir: 39 ************************************************************************************************* 3.2.2 MÉTODOS ITERATIVOS PARA A SOLUÇÃO DE SISTEMAS LINEARES No caso de sistemas lineares maiores e esparsos (com grande quantidade de zeros), os métodos diretos consomem grande memória computacional ou até podem não ser eficientes. Por isso, existem métodos de resolução iterativa de sistemas, ou seja, calculando uma solução inicial e outras sucessivamente, até que o critério de convergência seja satisfeito. Eles não serão aqui estudados, mas seguem links para o aluno que quiser se aprofundar no assunto. MÉTODO DE GAUSS-JACOBI (OU MÉTODO DE JACOBI OU MÉTODO DOS DESLOCAMENTOS SIMULTÂNEOS) http://vtp.ifsp.edu.br/nev/Sistema-jacobi/sistemajacobi.php? MÉTODO DE GAUSS-SEIDEL (OU MÉTODO DOS DESLOCAMENTOS SUCESSIVOS) http://vtp.ifsp.edu.br/nev/Sistema-seidel/sistemaseidel.php? 40 4 OPERAÇÕES COM MATRIZES 4.1 DETERMINANTE DE MATRIZ O cálculo de determinantes de matrizes pode não ter aplicação direta nos problemas de engenharia, mas aparece indiretamente para verificar se os sistemas lineares têm solução. Dado um sistema linear nxn e sua matriz de coeficientes A: - se det(A) ≠ 0 sistema tem solução única (possível e determinado) - se det(A) = 0 sistema tem infinitas soluções (possível e indeterminado) ou não tem soluções (impossível) Segundo Massago (2014), o determinante de uma matriz triangular (todos os elementos acima ou abaixo da diagonal são zeros) é dado pelo produto dos elementos das diagonais. O autor dá um exemplo: Sendo assim, o determinante de uma matriz qualquer pode ser calculado usando o escalonamento. Apartir do escalonamento de uma matriz A via MÉTODO DA ELIMINAÇÃO DE GAUSS da forma como foi estudado, determina-se uma matriz escalonada A‟ triangular superior equivalente. Para efeito do determinante, multiplicar uma linha por um múltipli de outra linha não causa alteração. Porém, trocar uma linha por outra inverte o sinal do determinante. Sendo assim, considerando que, na operação de escalonamento com pivoteamento parcial, trocam-se linhas “t” vezes: det [A] = (-1)t * (multiplicação dos elementos diagonais de A‟) Por exemplo, deseja-se calcular o determinante da seguinte matriz: A = ESCALONAMENTO COM PIVOTEAMENTO PARCIAL Como o pivô é 4, troca-se a linha 1 com a linha 3: 41 Realizando o escalonamento, chega-se a: Como o pivô é 0,25, não há troca de linhas. Realizando o escalonamento, chega-se a: Ao final do escalonamento, houve 1 troca de linhas (t=1). Sendo assim: det [A] = (-1)t vezes (multiplicação dos elementos diagonais de A‟) = (-1)1*(4*0,25*4) = -4 ************************************************************************************************* EXERCÍCIOS 4.1) Utilizando o esquema de entrada de dados do programa que calcula sistemas lineares por ELIMINAÇÃO DE GAUSS, crie um programa para calcular determinantes de matrizes. O resultado deve ficar como o do exemplo a seguir. 42 O algoritmo para multiplicar os elementos da diagonal de uma matriz pode ser: 1: det=1; 2: para j=1 até n { //para cada linha 3: det=det*a(j,j); 4: } OBS: lembre-se de contar as trocas de linha e multiplicar o determinante por (-1)t 4.2) No programa que calcula sistemas lineares por ELIMINAÇÃO DE GAUSS, incorpore uma verificação que pare a execução quando ele for possível e indeterminado ou impossível (det = 0). Teste o seguinte exemplo: ************************************************************************************************* 4.2 MULTIPLICAÇÃO DE MATRIZES Considere as matrizes a seguir e sua multiplicação: A = [auv]mxn B = [brs]nxp A*B = [ciJ]mxp Sendo assim: 43 (i=1,...,m; j=1,...,p) Um algoritmo para multiplicar matrizes é: 1: Para i=1 até m faça //para cada linha da matriz final (m é o nº de linhas da matriz A) 2: Para J=1 até p faça //para cada coluna da matriz final (p é o nº de colunas da matriz B) 3: c(i,J) = 0 4: Para k=1 até n faça // (n é o nº de linhas da matriz B) 5: c(i,J) = c(i,J)+a(i,k)*b(k,J) 6: Fim do laço 7: Fim do laço 8: Fim do laço 4.3 MATRIZ TRANSPOSTA Considere uma matriz A = a[i,J]mxn. Sua transposta será B = b[J,i]nxm. Um algoritmo para calcular a transposta de uma matriz é: 1: Para i=1 até m faça //para cada linha da matriz A 2: Para J=1 até n faça //para cada coluna da matriz A 3: b(J,i) = a(i,J) 4: Fim do laço 5: Fim do laço 4.4 MATRIZ INVERSA Um sistema linear pode ser resolvido numericamente por meio de um dos métodos apresentados no capítulo 3. Porém, poderia ser resolvido um sistema de equações usando a matriz inversa: [A]*[X]=[B] [X]=[A]-1*[B] Porém, é consenso que, computacionalmente, a obtenção da matriz inversa gera mais operações do que a resolução do sistema pelos outros métodos estudados. Isso será explícito no exemplo a seguir. Sendo assim, evita-se usar a matriz inversa em resoluções computacionais. Porém, se você se deparar com algum problema específico que exija a resolução de uma matriz inversa, Quadros e Bortoli (2009) apresentam uma solução, com um exemplo reproduzido seguir, semelhante ao processo de escalonamento de Gauss, porém usando a matriz aumentada com a matriz identidade (todas as diagonais valendo 1 e os outros elementos valendo zero). Ou seja, o que se faz é escalonar a matriz até que o sistema se torne triangular superior. Depois disso resolvem-se “n” sistemas triangulares superiores sendo que cada solução é uma coluna da matriz inversa (“n” é a ordem da matriz original). Resumindo, temos o seguinte esquema, segundo Castilho (2001): - matriz a inverter: 44 - matriz aumentada: - sistemas a resolver: - matriz inversa: Na prática, existem 2 maneiras de resolver o problema: 1) escalona-se a matriz aumentada e resolvem-se n sistemas triangulares, sendo o vetor “b” cada coluna da matriz identidade escalonada (atenção). 2) montar os sistemas sem escalonar, sendo “b” exatamente cada coluna da matriz identidade e escalonar cada sistema antes de resolver (como mostrado no esquema acima) Claramente o processo 1) resolve apenas 1 escalonamento e o 2) resolve “n” escalonamentos. Dependendo da ordem do sistema inicial, o processo 2) se torna inviável. Sendo assim, aqui, iremos resolver pelo processo 1). Então, da mesma maneira feita no processo de pivoteamento de Gauss, transforma-se a matriz aumentada em triangular superior antes de resolver cada sistema em separado. Castilho dá um exemplo numérico: matriz aumentada com I após escalonamento Resolve-se os 3 sistemas (triangulares superiores): 3 3 4 )1(X 45 E, portanto, a inversa fica: Castilho (2001) explica que, no exemplo anterior, a inversa de A é a própria matriz A, sendo que é um caso muito útil na verificação do método. Quadros e Bortoli (2005) apresentam outro exemplo: Resolvendo os sistemas triangulares superiores: Verificando os procedimentos anteriores, fica claro que o número de operações para encontrar a matriz inversa de um sistema é maior do que o custo para encontrar diretamente sua solução, pelo Método da Eliminação de Gauss. Isto sem considerar que, depois de encontrar a matriz inversa, deve-se multiplicá-la pelo vetor solução para definitivamente encontrar o vetor resposta do problema. Sendo assim, somente se usa o cálculo numérico de matriz inversa em casos onde é necessária sua dedução, o que não ocorre em sistemas lineares. ************************************************************************************************* 1 2 1 )2(X 5 6 6 )3(X 46 EXERCÍCIOS 4.3) Calcular, em papel, a inversa da matriz a seguir, pelo método apresentado anteriormente. 4.4) Partindo do algoritmo que calcula determinantes de matrizes, faça as adaptações e crie um programa para calcular inversas de matrizes. O resultado deve ficar como o do exemplo a seguir. O algoritmo para acrescentar a matriz identidade à matriz original lida pode ser: 1: para i=1 até n { //linhas 2: para j=n+1 até 2*n { //colunas 3: se i=(j-n) { //diagonal da identidade 4: a(i,j)=1; 5: } senão { 6: a(i,j)=0; 7: } fim do condicional 8: } fim do para i 9: } fim do para j ************************************************************************************************* 47 5 ZEROS DE FUNÇÕES Antes de iniciarmos os estudos deste capítulo, vamos relembrar as definições de FUNÇÃO e EQUAÇÃO. f(x) = x² + 2x – 8 FUNÇÃO: para qualquer valor de x, podemos calcular o valor correspondente de f(x). Especialmente quando xf =0, o valor de x é chamado raiz da função ou zero da função. Pode haver várias raízes, inclusive. x² –2x – 8 = 0 EQUAÇÃO: x só podem ser as raízes da função correspondente. x² –2x – 8 < 0 INEQUAÇÃO: existe um intervalo de valores de x que satisfazem a inequação. Sendo assim, este capítulo dedica-se ao estudo da obtenção de ZEROS DE FUNÇÕES, ou seja, raízes das EQUAÇÕES. Veja a figura que ilustra as raízes de uma função genérica: Figura: Raízes de uma função genérica. Fonte: Asano e Colli (2009). Sabemos que existe, por exemplo, a fórmula de Bhaskara para encontrar as raízes de polinômios de segundo grau. Porém, para polinômios de graus superiores, usam-se métodos do tipo iterativo para encontrar as soluções. Para polinômios de grau “n”, “n” raízes existem, podendo as mesmas serem reais ou complexas, diferentes ou não. Pode haver ainda um número infinito de soluções, como os casos de senos e cossenos. 5.1 MÉTODO DA BISSECÇÃO (OU DA DICOTOMIA) O Método da Dicotomia é um dos métodos numéricos mais simples usados para achar as raízes de uma função. O primeiro passo é isolar a raiz x* dentro de um intervalo onde a função seja monótona (ou crescente ou decrescente). Veja, na figura a seguir, o exemplo de um intervalo [a0,b0] de uma função f crescente, contendo a solução x* quando f=0: Figura: Intervalo [a0,b0] de uma função f crescente, contendo a solução x*. Fonte: Asano e Colli (2009). 48 Veja, no exemplo anterior, que f(a0)<0 e f(b0)>0 (seria ao contrário se a função fosse decrescente), o que leva a f(a0)*f(b0)<0. Em seguida, passamos a cercar a raiz com intervalos, cada intervalo com a metade do tamanho do intervalo anterior, até que se chegue a um intervalo tão pequeno que se possa assumir que o centro dele é a solução, com pouco erro. Asano e Colli (2009) apresentam um exemplo com a função f(x)=x3-20: 1) Escolhemos a0 = 2 porque f(a0) < 0 (2 3-20=-12<0) e b0 = 3 porque f(b0) > 0 (3 3- 20=7>0). Este é o chamado TEOREMA DE BOLZANO (sempre deve acontecer f(a)*f(b)<0, senão não cruza o eixo e não há raiz no intervalo). 2) Escolhemos o ponto médio do intervalo, ao qual chamaremos provisoriamente de c0: c0=2,5 3) Calculamos f(c0) = 2,5 3-20=-4,375<0. 4) Para descobrir se a raiz está à esquerda ou à direita de c0: - se a função for crescente (ou seja, f(a0) < 0) e f(c0) < 0 raiz à DIREITA de c0 - se a função for crescente (ou seja, f(a0) < 0) e f(c0) > 0 raiz à ESQUERDA de c0 - se a função for decrescente (ou seja, f(a0) > 0) e f(c0) < 0 raiz à ESQUERDA de c0 - se a função for decrescente (ou seja, f(a0) > 0) e f(c0) > 0 raiz à DIREITA de c0 Portanto, independentemente da função ser crescente ou decrescente, faça: - se f(a0)*f(c0) < 0 a raiz está à esquerda de c0 - se f(a0)*f(c0) > 0 a raiz está à direita de c0 No caso do exemplo, f(a0)*f(c0)=(-12)*(-4,375) > 0, portanto a raiz está à direita de c0. Assim, definimos um novo intervalo [a1,b1]= [c0,b0]= [2,5;3]. Portanto, a cada iteração “k”, deve-se verificar se a raiz está à esquerda ou à direita de ck. 4) Repete-se o procedimento até que o erro seja satisfatório. O erro real (ou seja, intervalo de confiança de x) en é dado por , assumindo que a solução está no ponto médio do intervalo. Outra forma de verificar a convergência é calcular f(c0), que deve ir se aproximando de zero a cada iteração. Abaixo está a tabela até a décima iteração, contendo a solução cn en. 49 Portanto, a raiz da função f(x)=x3-20 está no intervalo 2,71435 0,00049. Segundo Quadros e Bortoli (2009), este método possui baixa velocidade de convergência, mas a convergência é garantida. Não iremos explorar o algoritmo deste método ou realizar exercícios, visto ser mais interessante aprofundar o Método da Posição Falsa, uma derivação deste, com melhor velocidade de convergência. 5.2 MÉTODO DA POSIÇÃO FALSA Como comentado, este método é semelhante ao Método da Bissecção, só que mais elaborado. Ao invés de tomar a média aritmética entre “a” e “b”, o método da posição falsa toma a média aritmética ponderada entre “a” e “b” com pesos |f(b)| e |f(a)|, respectivamente. Então, considerando que f(a) e f(b) têm sinais opostos: Conforme Ruggiero (1996), graficamente, a solução x de cada iteração é dada por uma reta que passa pelos pontos (a,f(a)) e (b,f(b)) onde ela corta o eixo x: Fonte: Ruggiero (1996). 50 A sequência de operações, chamando x0=a, x1=b e x2=c neste caso, é: 1) Escolhemos x0 e x1 tais que f(x0)*f(x1)<0 (idêntico ao Método da Bissecção) 2) Aproximamos x2 a partir da expressão: 3) Se o critério de convergência não é satisfeito (aqui sendo f(x2) se aproximando de zero), segue-se. 4) Se f(x0)*f(x2)<0, mantém-se x0 e substitui-se x1 por x2 e retorne-se ao passo 2. Se f(x0)*f(x2)>0, mantém-se x1 e substitui-se x0 por x2 e retorne-se ao passo 2. (exatamente a mesma coisa que feito para o método da Dicotomia). Quadros e Bortoli (2009) apresentam um exemplo para a função p(x)=x3- 5x2+17x+21: ITERAÇÃO 0) 1) Escolhemos x0=-1 e x1=0 tais que f(x0)* f(x1)<0 (-2*21<0) 2) Aproximamos x2 a partir da expressão: 91304,0 23 21 21)2( 21*)1()2(*0 )0()1( )0(*)1()1(*0 )()( )(*)(* 10 1001 2 ff ff xfxf xfxxfx x 3) f(x2) = f(-0,91304) = (-0,91304) 3-5*(-0,91304)2+17*(-0,91304)+21=0,549. Veja que está bem longe de zero. Então parte para uma próxima iteração. Na tabela a seguir, encontram-se mais iterações até que haja uma convergência satisfatória: Obs: para cada iteração, manter a nomenclatura x0, x1 e x2. Portanto, uma das raízes da função é igual a -0,932115. Um algoritmo possível é, sendo dados a função f(x), o intervalo [a,b] tal que f(a)*f(b)<0, o erro alvo “ea” e o número máximo de iterações “itm”. Dado f(x), a, b, ea, itm 1: xa[0]=a; 2: xb[0]=b; 3: se ((f(xa[0])*f(xb[0]))>0) { //escolher outro intervalo 4: “Erro. Escolher outro intervalo com f(a)*f(b)<0” 51 5: } 6: k=0; //número da iteração 7: xm[k]=(xb[k]*f(xa[k])-xa[k]*f(xb[k]))/(f(xa[k])-f(xb[k])); //solução da iteração 0 8: se ((f(xa[k])*f(xm[k]))<=0) { //raiz da próxima iteração está à esquerda de xm 9: xa[k+1]=xa[k]; 10: xb[k+1]=xm[k]; 11: } senão { //raiz da próxima iteração está à direita de xm 12: xa[k+1]=xm[k]; 13: xb[k+1]=xb[k]; 14: } fim do condicional 15: En[k]=módulo(f(xm[k])); //erro relativo da iteração 0 16: enquanto (En[k]>ea) { //a partir da iteração 1 17: k=k+1; 18: se (k=itm) { 19: break; //sai do laço enquanto 20: } fim do condicional 21: xm[k]=(xb[k]*f(xa[k])-xa[k]*f(xb[k]))/(f(xa[k])-f(xb[k])); //solução da iteração k 22: se ((f(xa[k])* f(xm[k]))<=0) { //raiz da próxima iteração está à esquerda de xm 23: xa[k+1]=xa[k]; 24: xb[k+1]=xm[k]; 25: } senão { //raiz da próxima iteração está à direita de xm 26: xa[k+1]=xm[k]; 27: xb[k+1]=xb[k]; 28: } fim do condicional 29: En[k]=módulo(f(xm[k])); //erro relativo da iteração k 30: } fim do loop equanto O professor irá enviar aos grupos, via e-mail, um programa para calcular a raiz de uma função qualquer pelo Método da POSIÇÃO FALSA, contendo a plotagem do gráfico da função. A plotagem de gráficos necessita de arquivos específicos que estão na pasta “funcoes_grafico”. Tais arquivos não são necessários de serem estudados. Eles são funções disponibilizadas on-line, mas que não são de código aberto (não se pode editar), apesar de serem gratuitas. O aluno deve estudar somente os arquivos html e js, além de testar e entender o funcionamento do programa. Segue um exemplo de resultado: 52 ************************************************************************************************* EXERCÍCIO5.1) Calcular, em papel, a raiz da função no intervalo , usando o Método da Posição Falsa. Faça iterações até que a convergência seja menor que 0,001. 5.2) Calcular o exercício anterior usando o algoritmo passo-a-passo. 5.3) Considere uma viga biapoiada submetida a uma ação q que segue uma variação triangular conforme a figura a seguir. Neste caso, o diagrama de momento fletor é dado pela função M(x) e o diagrama de esforço cortante é dado por Q(x) (veja a seguir). Usando o Método da Posição Falsa, calcule a posição x onde o momento fletor é máximo (Q=0), para o caso de q = 1,5 tf/m e L = 4 m. Calcule o valor do momento máximo. Faça iterações até que o erro seja menor que 0,001. ************************************************************************************************* 53 5.3 MÉTODO DE NEWTON-RAPSON (OU DAS TANGENTES) A inspiração do método é vista na figura a seguir, onde olhamos para a reta tangente ao gráfico de f no ponto (xk,f(xk)) e definimos xk+1 como sendo o ponto de encontro dessa reta com o eixo x. Assim, sempre se terá uma aproximação mais perto da solução x*. Figura: Visualização geométrica do Método de Newton-Rapson. Fonte: Asano e Colli (2009). Deduzindo a expressão, a partir da figura anterior: )( )( )(' 1 kk k k xx xf xf (a derivada é o coeficiente angular da reta tangente) )()(*)(' 1 kkkk xfxxxf )(*)('*)(' 1 kkkkk xfxxfxxf )(' )(*)(' 1 k kkk k xf xfxxf x )(' )( 1 k k kk xf xf xx Segundo Quadros e Bortoli (2009), este método pode ser mais eficiente por causa da velocidade de convergência e da precisão. Porém, enquanto os métodos anteriormente apresentados só necessitam da resolução da função em si, este método necessita da resolução de sua derivada. Agora não é mais necessário definir um intervalo inicial, mas somente um valor de x inicial. Portanto, obviamente, não há verificação inicial. Asano e Colli (2009) apresentam o mesmo exemplo mostrado em capítulo anterior resolvido pelo Método da Posição Falsa, para encontrar uma raiz da função f(x) = x3-20: Chuta-se x0 = 3: 54 Observe, pela tabela anterior, que a convergência (f(xn) se aproximando de 0) foi rápida, encontrando-se a solução (x = 2,7144) em apenas 3 iterações, sendo que, pelo Método da Posição Falsa, precisou-se de 4 iterações (considerando a mesma convergência alvo = 0,001). Uma desvantagem do Método de Newon-Rapson é a necessidade da derivada da função. Então, como uma evolução, usa-se, ao invés da derivada direta, a derivação numérica, dando origem ao Método das Secantes. ************************************************************************************************* 5.4 MÉTODO DAS SECANTES Como comentado, a partir do Método de Newton-Rapson, no lugar da derivada direta da função, usa-se a derivação numérica, dada por: Sendo assim, o termo geral de iterações é dado por: )(' )( 1 k k kk xf xf xx Segundo Quadros e Bortoli (2009), este método é quase tão eficiente quanto o Método de Newton-Rapson. Os autores apresentam o exemplo a seguir, que resolve o polinômio utilizando . Atente-se que você deve começar os cálculos a partir da iteração i=1 e entrar com x0 e x1 (DUAS APROXIMAÇÕES, em que xi é a solução), porque ela já vai usar xi-1 = x0 na primeira iteração. Considera-se o teste de convergência como o próprio f(xi), que deve se aproximar do zero na convergência para a raiz. ************************************************************************************************* EXERCÍCIOS 5.4) Calcular, em papel, a raiz da função usando o Método das Secantes (a mesma função do exercício 5.1). Usar x0=0 e x1=1. Faça iterações até que a convergência seja menor que 0,001. 55 5.5) Adapte o algoritmo do Método da Posição Falsa para que resolva a raiz de uma função pelo Método das Secantes. 5.6) Duplique o programa do Método da Posição Falsa e adapte-o para que resolva a raiz de uma função pelo Método das Secantes. A aparência deve ser semelhante, mas com x0 e x1 (aproximações iniciais). Veja: ************************************************************************************************* 56 Quadros e Bortoli (2009) citam que os algoritmos de resolução de zeros de funções devem prever a possibilidade de ocorrência de divisão por zero ou por um valor muito pequeno, o que pode baixar muito a velocidade de convergência. Citam ainda os autores que pode haver comportamento oscilante nas iterações quando a função não apresenta raízes reais, quando a estimativa inicial não for adequada ou quando a função for simétrica em relação ao eixo x (f‟‟(x*)=0). x* é a raiz da função. Quadros e Bortoli (2009) apresentam em seu trabalho mais alguns métodos para encontrar zeros de funções, que podem ser usados em estudos mais aprofundados. 57 6 FUNÇÕES DE INTERPOLAÇÃO DE DADOS Frequentemente nos deparamos com um conjunto discreto de valores que podem ser dados na forma de tabela ou de um conjunto de medidas. Por exemplo, temos o seguinte trecho de tabela dada pela norma de cálculo de ação de vento em edificações, para a obtenção de um fator chamado S2: Figura: Trecho de tabela para a definição do fator S2. Fonte: NBR 6123 (1988). Veja, na tabela anterior, que, se tivermos um valor de z igual a 22, não temos como encontrar diretamente S2, seja qual for a classe e a categoria. Porém, existe S2 para z = 20 e para z = 30. E o 22 fica entre esses 2 valores. Neste caso, temos que fazer uma interpolação. Fazer a conta da interpolação para um único valor é fácil. Porém, nos casos práticos, temos muitos valores a interpolar, o que acaba necessitando interpolar todos os pontos e encontrar uma função contínua que satisfaça todos os pontos. É a chamada função interpoladora. Figura: Polinômio interpolador para alguns dados discretos. Fonte: http://www.math.ist.utl.pt/~calves/cursos/Interpola.HTM. Acesso em 04/11/2014. A função interpoladora pode ser polinomial, exponencial, entre outras. Observe que a função interpoladora somente deve ser usada dentro do domínio de pontos conhecidos (x0 a x3 no caso da figura anterior). O domínio fora dos pontos conhecidos deve ser tratado como uma extrapolação, que será estudada em capítulo posterior. 58 6.1 INTERPOLAÇÃO POLINOMIAL Considere que se deseja interpolar n+1 pontos (x,f(x)) por um polinômio p(x), de grau menor ou igual a n. Sendo assim: O sistema apresentado tem n+1 equações e incógnitas (a0 até an). Para resolver o sistema, é necessário que a matriz dos coeficientes (a seguir) tenha seu determinante diferente de zero para que a solução seja única e determinada. Matricialmente, o sistema apresentado fica: )( )( )( * 1 1 1 1 00 1 0 2 2 11 2 00 nn a n n n n nn xf xf xf a a a x x x xx xx xx Computacionalmente, pode-se resolver o sistema por algum dos métodos explicados em capítulo anterior. Quadros e Bortoli (2009) apresentam um exemplo, em que se pede para encontrar um polinômio de grau 2 para interpolar os seguintes pontos: Adotando , tem-se: p2(x0) = f(x0) a0+a1*(-2)+a2*(-2)
Compartilhar