Universidade Federal de Uberl^ andia Faculdade de Matem atica Cerico alculo Num Prof. Jos e Eduardo Castilho Marco de 2001 Conte´udo 1 Introdu¸c˜ao 1 1.1 OMatLab .................................... 3 1.1.1 C´alculonaJaneladeComandos ..................... 3 1.1.2 M-arquivos ................................ 7 1.2 Exerc´icios ..................................... 10 2 Zeros de Fun¸c˜oes 11 2.1 IsolamentodasRa´izes .............................. 11 2.2 Refinamento .................................... 14 2.3 M´etododaBissec¸c˜ao ............................... 14 2.3.1 EstudodaConvergˆencia ......................... 15 2.3.2 EstimativadoN´umerodeItera¸c˜oes . . . . . . . . . . . . . . . . . . . 16 2.4 M´etodoIterativoLinear(M.I.L.) ........................ 17 2.4.1 Crit´eriodeParada ............................ 19 2.5 M´etododeNewton-Raphson(M.N.R)...................... 20 2.6 OrdemdeConvergˆencia ............................. 22 2.7 Observa¸c˜oesFinais ................................ 24 2.8 Exerc´icios ..................................... 25 3 Sistemas Lineares 27 3.1 M´etodosDiretos.................................. 27 3.1.1 SistemaTriangularSuperior ....................... 28 3.1.2 M´etododeElimina¸c˜aodeGauss ..................... 29 3.1.3 PivotamentoParcial ........................... 32 3.1.4 C´alculodaMatrizInversa ........................ 34 3.2 M´etodosIterativos ................................ 38 3.2.1 Crit´eriodeConvergˆencia......................... 40 3.2.2 M´etodoIterativodeGauss-Jacobi . . . . . . . . . . . . . . . . . . . . 40 3.2.3 Crit´eriodasLinhas ............................ 42 3.2.4 M´etodoIterativodeGauss-Seidel . . . . . . . . . . . . . . . . . . . . 44 3.2.5 Crit´eriodeSassenfeld........................... 45 3.3 Observa¸c˜oesFinais ................................ 46 i ´ CONTEUDO ii 3.4 Exerc´icios ..................................... 47 4 Ajuste de Curvas: M´etodo dos M´inimos Quadrados 49 4.1 M´etododosM´inimosQuadrados -CasoDiscreto. . . . . . . . . . . . . . . . 49 4.2 M´etodo dos M´inimos Quadrados -Caso Cont´inuo . . . . . . . . . . . . . . . 53 4.3 AjusteN˜aoLinear ................................ 55 4.4 Observa¸c˜oesFinais ................................ 56 4.5 Exerc´icios ..................................... 58 5 Interpola¸c˜ao Polinomial 60 5.1 FormadeLagrange ................................ 62 5.2 FormadeNewton ................................. 63 5.2.1 Constru¸c˜aodoPolinˆomio......................... 64 5.3 EstudodoErro .................................. 65 5.4 EscolhadosPontos ................................ 67 5.5 Interpola¸c˜aoInversa ............................... 67 5.6 Observa¸c˜oesFinais ................................ 68 5.7 Exerc´icios ..................................... 70 6 Integra¸c˜ao Num´erica -F´ormulas de Newton Cˆotes 72 6.1 RegradoTrap´ezio ................................ 72 6.2 C´alculodoErro ................................. 73 6.3 RegradoTrap´ezioRepetida ........................... 75 6.4 RegradeSimpson1/3 .............................. 76 6.5 RegradeSimpsonRepetida ........................... 78 6.6 Observa¸c˜oesFinais ................................ 79 6.7 Exerc´icios ..................................... 79 7 Equa¸c˜oes Diferenciais Ordin´arias (P.V.I.) 81 7.1 M´etodoEuler ................................... 81 7.2 M´etodosdaS´eriedeTaylor ........................... 83 7.3 M´etodosdeRunge-Kutta............................. 85 7.4 M´etodosdeAdans-Bashforth .......................... 88 7.4.1 M´etodosExpl´icitos............................ 89 7.4.2 M´etodosImpl´icitos ............................ 90 7.5 Equa¸c˜oesdeOrdemSuperior........................... 92 7.6 Exerc´icios ..................................... 93 Cap´itulo 1 Introdu¸c˜ao O C´alculo Num´erico tem por objetivo estudar esquemas num´ericos (algoritmos num´ericos) para resolu¸c˜ao de problemas que podem ser representados por um modelo matem´atico. Um esquema ´e eficiente quando este apresenta solu¸c˜oes dentro de uma precis˜ao desejada com custo computacional (tempo de execu¸c˜ao + mem´oria) baixo. Os esquemas num´ericos nos fornecem aproxima¸c˜oes para o que seria a solu¸c˜ao exata do problema. Os erros cometidos nesta aproxima¸c˜ao s˜ao decorrentes da discretiza¸c˜ao do problema, ou seja passar do modelo matem´atico para o esquema num´erico, e da forma como as m´aquinas representam os dados num´ericos. Como exemplo de discretiza¸c˜ao consideremos que desejamos calcular uma aproxima¸c˜ao para a derivada de uma fun¸c˜ao f(x) num ponto ¯x. O modelo matem´atico ´e dado por f(¯x) = lim f(¯x + h) - f(¯x) h0 h ! Um esquema num´erico para aproximar a derivada ´e dado por tomar h “pequeno” e calcular f0(¯x) ˜ f(¯x + h) - f(¯x) h Neste caso quanto menor for o valor de h mais preciso ser´a o resultado, mas em geral, este esquema n˜ao fornecer´a a solu¸c˜ao exata. A representa¸c˜ao de n´umeros em m´aquinas digitais (calculadoras, computadores, etc) ´e feita na forma de ponto flutuante com um n´umero finito de d´igito. Logo os n´umeros que tem representa¸c˜ao infinita (Ex. 1=3, ¼, p2) s˜ao representados de forma truncada. Com isto algumas das propriedades da aritm´etica real n˜ao valem na aritm´etica computacional. Como exemplo, na aritm´etica computacional temos nn. ak = 1 . ak; NN k=0 6k=0 onde estamos considerando que no primeiro somat´orio para cada k fazemos ak=N e depois somamos e no segundo somat´orio somamos todos os ak e o resultado da soma dividimos por 1 ˜ CAP´ITULO 1. INTRODUC¸ AO 2 N. Do ponto de vista anal´itico, as duas express˜oes s˜ao equivalentes, mas a segunda forma apresenta melhor resultado do ponto de vista computacional, pois realiza menos opera¸c˜oes e comete menos erro de truncamento. Outro exemplo interessante ´e que em aritm´etica computacional ´e poss´ivel que para um dado A exista um e = 0 tal que A + e = A. Analiticamente a express˜ao acima ´e verdadeira se e somente se e = 0. Outro fator que pode influenciar no resultado ´e o tipo de m´aquina em que estamos trabalhando. Numa calculadora simples que represente os n´umeros com 7 d´igito ter´iamos 1=3+1=3+1=3=0:9999999 Enquanto que calculadoras mais avan¸cadas ter´iamos como resposta um falso 1, pois na realidade, internamente estas calculadoras trabalham com mais d´igito do que ´e apresentado no visor e antes do resultado ser apresentado este ´e arredondado. Os esquemas num´ericos s˜ao classificados como esquemas diretos e esquemas iterativos. Os esquemas diretos s˜ao aqueles que fornecem a solu¸c˜ao ap´os um n´umero finito de passos. Por exemplo o esquema que apresentamos para o c´alculo da derivada. Os esquemas iterativos s˜ao aqueles que repetem um n´umero de passos at´e que um crit´erio de parada seja satisfeito. Como exemplo considere o algoritmo que ´e usado para determinar a precis˜ao de uma m´aquina digital Algoritmo: Epsilon da M´aquina Ep 1 . Enquanto (1 + Ep) > 1, fa¸ca: Ep Ep=2 . fim enquanto OutPut: 2Ep O crit´erio de parada ´e (1 + Ep) = 1 e cada execu¸c˜ao do la¸co Enquanto ´e chamado de itera¸c˜ao. M´aquinas diferentes apresentar˜ao resultados diferentes. Um outro fator que pode influenciar nos resultados ´e a linguagem de programa¸c˜ao usada na implementa¸c˜ao dos algoritmos (Pascal, Fortran, C++, MatLab , etc). Diferentes linguagens podem apresentar diferentes resultados. E mesmo quando usamos uma mesma linguagem, mas compiladores diferentes (Ex. C++ da Borland e C++ da Microsoft), os resultados podem apresentar diferen¸cas. Existem v´arias bibliotecas de rotinas num´ericas em diversas linguagens e algumas dispon´iveis na Internet. Um exemplo ´e a LIMPACK: uma cole¸c˜ao de rotinas em Fortran para solu¸c˜ao de sistemas lineares. Para exemplificar os esquemas num´ericos, que estudaremos nos pr´oximos cap´itulo, usaremos o MatLab pela sua facilidade de programa¸c˜ao. Todo o material desta apostila ´e baseado nas referencias: [?]e[?]. ˜ CAP´ITULO 1. INTRODUC¸ AO 3 1.1 O MatLab O MatLab surgiu nos anos 1970 como um Laborat´orio de Matrizes para auxiliar os cursos ´ de Teoria Matricial, Algebra Linear e Analise Num´erica. Hoje, a capacidade do MatLab se estende muito al´em da manipula¸c˜ao de matrizes. Ele ´e tanto um ambiente quanto uma linguagem de programa¸c˜ao, e um de seus aspectos mais poderosos ´e que os problemas e as solu¸c˜oes s˜ao expressos numa linguagem matem´atica bem familiar. Devido a sua capacidade de fazer c´alculos, visualiza¸c˜ao gr´afica e programa¸c˜ao, num ambiente de f´acil uso, o MatLab torna-se uma ferramenta eficiente para a compreens˜ao tanto de t´opicos fundamentais quanto avan¸cados a uma gama de disciplinas. Nosso objetivo ´e dar uma r´apida vis˜ao dos comandos e fun¸c˜oes b´asicas do MatLab para exemplificar os t´opicos do curso de C´alculo Num´erico. Maiores detalhes podem ser obtidos em ??. Apesar das ultimas vers˜oes do MatLab ter expandido sua capacidade, o elemento b´asico dos dados ainda ´e um vetor, o qual n˜ao requer declara¸c˜ao de dimens˜ao ou tipo de vari´avel. O MatLab ´e um sistema interativo, onde os comandos podem ser executados na janela de comandos ou por programas. Estes programas s˜ao conhecidos como m-arquivos ( ou arquivos com extens˜ao .m) e ser˜ao discutidos posteriormente. Em primeiro lugar vamos discutir alguns comandos b´asicos que ser˜´c˜ ao util para a manipula¸ao de dados na janela de comandos e nos m-arquivos. 1.1.1 C´alculo na Janela de Comandos Um c´alculo simples pode ser executado na janela de comandos digitando as instru¸c˜oes no prompt como vocˆe faria numa calculadora. Por exemplo EDU>> 3*4 +5 ans = 17 o resultado ´e mostrado na tela como ans ( abreviatura de “answer”). Os s´imbolos dos operadores aritm´eticos s˜ao dados na Tabela 1.1. As express˜oes s˜ao calculadas da esquerda para a direita, com a potencia¸c˜ao tendo a maior precedˆencia, seguido da multiplica¸c˜ao e divis˜ao (mesma precedˆencia) e pela adi¸c˜ao e subtra¸c˜ao (tamb´em com mesma precedˆencia). As Vari´aveis A forma de armazenar o resultado para uso posterior ´e pelo uso de vari´aveis. EDU>> s=3+4+7+12 s= ˜ CAP´ITULO 1. INTRODUC¸ AO 4 Tabela 1.1: Operadores Aritm´eticos Opera¸c˜ao S´imbolo Adi¸c˜ao a + b Multiplica¸c˜ao a * b Subtra¸c˜ao a - b Divis˜ao a=b ou bna Potencia¸c˜ao abb O nome da vari´avel pode consistir de no m´aximo 31 caracteres, iniciando sempre por um caracter alfa seguido de qualquer combina¸c˜ao de caracteres do tipo alfa , num´erico e underscores. Ex. resultado_da_soma_2. Ao contr´ario de outras linguagens, o MatLab diferencia as vari´aveis que usam letras min´usculas e mai´usculas. Isto ´e as vari´aveis Contas, contas, conTas e CoNtAs, s˜ao consideradas como quatro vari´aveis diferentes. Todas as vari´aveis s˜ao armazenadas internamente e podem ser usadas a qualquer momento. Para saber quais as vari´aveis que est˜ao ativas utilizamos o comando who. Para eliminar a vari´avel conta, usamos o comando clear conta. As vari´aveis s˜ao tratadas como matrizes, apesar dos escalares n˜ao serem apresentados na nota¸c˜ao matricial. Um vetor linha pode ser definido como EDU>> x=[1 2 3 4] x= 1234 Tamb´em podemos separar os elementos por v´irgula. J´a um vetor coluna pode ser definido da seguinte forma EDU>> y=[5; 6; 7; 8] y= 5 6 7 8 Um exemplo de uma matriz 3 × 4. EDU>> a=[1234;5678;9101112] a = 1 2 3 4 5 6 7 8 9 10 11 12 ˜ CAP´ITULO 1. INTRODUC¸ AO 5 Os elementos na i-´esima linha e na j-´esima coluna, denotados por aij podem ser obtidos pelo comando a(i,j), por exemplo a(2,3)=7. Note que os elementos no MatLab s˜ao indexados iniciando em 1. Em algumas situa¸c˜oes necessitamos de vetores com alguma estrutura particular. Por exemplo, um vetor cujo o primeiro termo vale ¡2 e o ultimo vale 3 e os termos intermedi´arios variam um passo de 0:5. Este vetor pode ser definido pela linha de comando EDU>> v=-2:0.5:3 v= Columns 1 through 7 -2.0000 -1.5000 -1.0000 -0.5000 0 0.5000 1.0000 Columns 8 through 11 1.5000 2.0000 2.5000 3.0000 De uma forma geral este comando tem a sintaxe v=a:passo:b. Quando o passo ´e igual a um podemos escrever o comando na forma reduzida v=a:b. Algumas matrizes elementares podem ser geradas atrav´es de comandos simples, por exemplo: A matriz identidade: EDU>>I=eye(3) I= 100 010 001 Matriz n × m formada por 1’s: EDU>> A=ones(2,3) A= 111 111 Matriz Nula de ordem n × m: EDU>> B=zeros(3,4) B= 0000 0000 0000 ˜ CAP´ITULO 1. INTRODUC¸ AO 6 Opera¸c˜oes com Matrizes e Vetores As opera¸c˜oes de subtra¸c˜ao,adi¸c˜ao e multiplica¸c˜ao entre matrizes s˜ao definidas como o usual. O s´imbolo da divis˜ao representa o c´alculo da multiplica¸c˜ao pela inversa da matriz, por exemplo EDU>> A=[1 2 1;3 2 4;5 3 2]; EDU>> A/A ans = 1 0 0 0 1 0 0 0 1 Note que neste exemplo terminamos o comando que define a matriz A com um ponto e v´irgula. Isto faz com que o resultado do comando (ou de uma express˜ao em geral) n˜ao seja apresentado no monitor. Isto e util para programas de grande porte, pois este processo ´´ de apresentar os resultados no monitor consome muito tempo de execu¸c˜ao. Entre vetores e matrizes de mesma dimens˜ao ´e poss´ivel operar elemento com elemento. Isto ´e poss´ivel atrav´es de uma nota¸c˜ao especial, que consiste de usar um ponto antes do s´imbolo da opera¸c˜ao. Na Tabela 1.2 damos um resumo destas opera¸c˜oes aplicada a dois vetores a e b de dimens˜ao n. Tabela 1.2: Opera¸c˜oes Elementares entre Vetores Opera¸c˜ao S´imbolo Resultado Adi¸c˜ao a+b [a1 + b1, a2 + b2, a3 + b3, . . . , an + bn] Subtra¸c˜ao a-b [a1 - b1, a2 - b2, a3 - b3, . . . , an - bn] Multiplica¸c˜ao a.*b [a1 * b1, a2 * b2, a3 * b3, . . . , an * bn] Divis˜ao a./b [a1=b1, a2=b2, a3=b3, . . . , an=bn] Potencia¸c˜ao a:bb [a1 b1 , a2 b2 , a3 b3 , . . . , an bn ] Como na maioria das linguagens de programa¸c˜ao, o MatLab oferece diversas fun¸c˜oes elementares que s˜ao importantes em matem´atica. A Tabela 1.3 apresenta uma lista destas fun¸c˜oes e sua sintaxe. Gr´aficos Para plotar um gr´afico no MatLab , devemos criar dois vetores de mesma dimens˜ao x e f, onde x corresponde aos valores do eixo x e f os valores da fun¸c˜ao nestes pontos. O ˜ CAP´ITULO 1. INTRODUC¸ AO 7 Tabela 1.3: Fun¸c˜oes Elementares Fun¸c˜ao Sintaxe Valor Absoluto abs(x) Arco Co-seno acos(x) Arco Seno asin(x) Co-seno cos(x) Exponencial ex exp(x) Logaritmo Natural log(x) Logaritmo base 10 log10(x) Seno sin(x) Raiz Quadrada sqrt(x) Tangente tan(x) gr´afico ´e gerado pelo comando plot(x,f). Se desejamos gerar o gr´afico da fun¸c˜ao sen(x) no intervalo [¡¼, ¼] devemos proceder da seguinte forma: EDU>> x=-pi:0.01:pi; EDU>> f=sin(x); EDU>> plot(x,f) Note, que na defini¸c˜ao do vetor x, usamos o passo igual a 0:01. Isto determina a quantidade de pontos que o comando plot usa par gerar o gr´afico. Quanto mais pontos mais perfeito ser´a o gr´afico (em contra partida maior o tempo de execu¸c˜ao). se tiv´essemos usado o passo 0:5 n˜ao ter´iamos um gr´afico de boa qualidade. 1.1.2 M-arquivos Existem dois tipos de programas em MatLab : scripts e funtions. Ambos devem ser salvos com extens˜ao .m no diret´orio corrente. Uma diferen¸ca b´asica entre os dois ´e que os scripts trata as vari´aveis, nele definidas, como vari´aveis globais, enquanto as functions trata as vari´aveis como vari´aveis locais. Desta forma a functions tem que ter um valor de retorno. Scripts Os scripts permite que um conjunto de comandos e defini¸c˜oes sejam executados atrav´es de um ´ unico comando na janela de comandos. Como exemplo, o script seguinte calcula a aproxima¸c˜ao da derivada de sen(x) usando diferen¸cas finitas. % Aproximacao da derivada do seno % Usando o operardor de diferen\c{c}a finita progressiva. ˜ CAP´ITULO 1. INTRODUC¸ AO 8 clear; h=0.0001; x=input('Entre com o valor de, x='); % Atribui Valores a x disp('O valor da aproximacao eh...') % Mostra mensagem no monitor dsen=(sin(x+h)-sin(x))/h As primeiras duas linha s˜ao coment´arios que descrevem o script. Na quinta linha temos o comando que permite atribuir valores a uma vari´avel. E na sexta linha o comando que permite mostrar uma mensagem no monitor. Vamos supor que este arquivo seja salvo com o nome de devira_seno.m. Para executar o script digitamos seu nome no prompt do MatLab . Functions Numa fun¸c˜ao em MatLab a primeira linha ´e da forma function y=nome(argumentos). A fun¸c˜ao se troca informa¸c˜oes com o MatLab workspace por interm´edio da vari´avel y e dos argumentos. Para ilustrar o uso de fun¸c˜oes em MatLab considere o seguinte c´odigo function dsen=deriva_seno(x,h) % Aproximacao da derivada do seno % Usando o operardor de diferen\c{c}a finita progressiva. dsen=(sin(x+h)-sin(x))/h; Apesar deste arquivo poder ser salvo com um nome qualquer, ´e usual usar o mesmo nome da fun¸c˜ao, ou seja, deriva_seno.m. Para execut´a-lo devemos digitar seu nome e informar os valores dos argumentos, por exemplo, EDU>>y= deriva_seno(3.14,0.001) o que forneceria em y uma aproxima¸c˜ao da derivada da fun¸c˜ao seno em 3:14 Uma diferen¸ca importante entre esta vers˜ao, usando function, com a anterior ´e que o valor calculado pode ser atribu´ido a uma vari´avel. Al´em disso, agora podemos escolher o valor de h, que na vers˜ao anterior estava fixo em h=0.0001. Vale notar que no primeiro caso todas as vari´aveis do scrips est˜ao ativas, isto ´e s˜ao vari´aveis globais. Enquanto que no segundo ca so as vari´aveis s˜ao locais, isto ´e, a vari´avel h s´o ´e ativa na execu¸c˜ao da fun¸c˜ao Controle de Fluxo O controle de fluxo ´e um recurso que permite que resultados anteriores influenciem opera¸c˜oes futuras. Como em outras linguagens, o MatLab possui recursos que permitem o controle de fluxo de execu¸c˜ao de comandos, com base em estruturas de tomada de decis˜oes. Apresentamos as estrutura de loops for, loops while e if-else-end. A forma geral do loop for ´e ˜ CAP´ITULO 1. INTRODUC¸ AO 9 for x = vetor comandos... end Os comandos entre for e end s˜ao executados uma vez para cada coluna de vetor. A cada itera¸c˜ao atribui-se a x a pr´oxima coluna de vetor. Por exemplo EDU>> for n=1:5 x(n) = cos(n*pi/2); end EDU>> x x= 0.0000 -1.0000 -0.0000 1.0000 0.0000 Traduzindo, isto diz que para n igual a 1 at´e 10 calcule os comandos at´e end. Ao contr´ario do loop for, que executa um grupo de comandos um n´umero fixo de vezes, o loop while executa um grupo um de comandos quantas vezes forem necess´arias para que uma condi¸c˜ao seja negada. Sua forma geral ´e while expressao comandos... end O grupo de comandos entre while e end s˜ao executados at´e que a express˜ao assuma um valor falso. Por exemplo, EDU>> while abs(x(n)-x(n-1)) > 10^(-6) x(n) = 2*x(n-1) + 1/4; n=n+1; end Neste caso o grupo de comandos s˜ao executados at´e que o valor absoluto da diferen¸ca entre dois valores consecutivos seja menor ou igual a 10¡6 . A estrutura if-else-end permite que grupos de comandos sejam executados por um teste relacional. A forma geral ´e dada por if expressao comandos 1... else comandos 2... end Se a express˜ao for verdadeira ´e executado o grupo de comandos 1, caso contr´ario ´e executado o grupo de comandos 2. Esta estrutura permite o uso da forma mais simples que envolve s´o um condicional ˜ CAP´ITULO 1. INTRODUC¸ AO 10 if expressao comandos ... end Como exemplo considere o seguinte fragmento de c´odigo que calcula o valor absoluto de um n´umero if x < 0 x=-x; end Isto ´e, se x for menor que zero ent˜ao troca de sinal, caso contr´ario nada ´e feito. 1.2 Exerc´icios Exerc´icio 1.1 Usando o esquema num´erico para a aproxima¸c˜ao da derivada dado abaixo ache uma aproxima¸c˜ao para f0(¼), onde f(x)= sen(x) e tome h =0:1, 0:01, 0:001;::. 10¡10 . Repita os c´alculos para f0(0). Comente os resultados. f0(¯x) ˜ f(¯x + h) - f(¯x) h Exerc´icio 1.2 Fa¸ca um programa que calcule 107 A + . 10¡7 k=1 com A = 10, 102 , 103 ;:::, 1015 . Comente os resultados. Exerc´icio 1.3 Calcule a precis˜ao de sua m´aquina usando o algoritmo Algoritmo: Epsilon da M´aquina Input: A : n´umero que represente a grandeza Ep 1 . Enquanto (A + Ep) > 1, fa¸ca: Ep Ep=2 . fim enquanto Output: Imprimir 2Ep tomando A =1, 10, 100, 1000. Comente os resultados. Cap´itulo 2 Zeros de Fun¸c˜oes Neste cap´itulo estudaremos esquemas num´ericos para resolver equa¸c˜oes da forma f(x) = 0. Na maioria dos casos estas equa¸c˜oes n˜ao tem solu¸c˜ao alg´ebrica como h´a para as equa¸c˜oes de 2 o ¯ grau. No entanto esquemas num´ericos podem fornecer uma solu¸c˜ao satisfat´oria. O processo para encontrar uma solu¸c˜ao envolve duas fases: Fase I Isolamento das ra´izes -Consiste em achar um intervalo fechado [a, b] que cont´em a raiz. Fase II Refinamento -Partindo de uma aproxima¸c˜ao inicial refinamos a solu¸c˜ao at´e que certos crit´erios sejam satisfeitos. 2.1 Isolamento das Ra´izes Um n´umero x que satisfaz a equa¸c˜ao f(x) = 0 ´e chamado de raiz ou zero de f. O objetivo ´e encontrar um intervalo [a, b], de pequena amplitude ( b - a . 1), que contenha a raiz que desejamos encontrar. Para isto usaremos duas estrat´egias: An´alise Gr´afica e Tabelamento da fun¸c˜ao. A an´alise gr´afica ´e baseada na id´eia de que, a partir da equa¸c˜ao f(x) = 0, podemos obter uma equa¸c˜ao equivalente g(x) - h(x) = 0, onde g e h sejam fun¸c˜oes mais simples e de f´acil an´alise gr´afica. Esbo¸cando o gr´afico de g e h podemos determinar os pontos x, onde as curvas se interceptam, pois estes pontos ser˜ao as ra´izes de f(x)( g(»)= h(») f(»)=0 ). . Exemplo 2.1.1 Sendo f(x)= e¡x - x temos f(x)= g(x) - h(x), onde g(x)= e¡x e h(x)= x. Na Figura 2.1 temos que as curvas se interceptam no intervalo [0, 1]. Tamb´em podemos observar que pelo comportamento das fun¸c˜oes g(x) e h(x) estas fun¸c˜oes n˜ao v˜ao se interceptar em nenhum outro ponto. Logo f(x) admite uma ´unica raiz. Na pr´atica usamos algum software matem´atico para esbo¸car os gr´aficos. Quanto menor for a amplitude do intervalo que cont´em a raiz, mais eficiente ser´a a Fase de Refinamento. Para obtermos um intervalo de menor amplitude usaremos a estrat´egia do tabelamento que ´e baseada no seguinte Teorema. 11 ˜ CAP´ITULO 2. ZEROS DE FUNC¸OES 12 Figura 2.1: Gr´aficos de g(x)e h(x) Teorema 2.1.1 Seja f(x) uma fun¸c˜ao cont´inua num intervalo [a, b]. Se f(a)f(b) < 0 ent˜ao existe pelo menos uma raiz . . [a, b]. O Teorema garante a existˆencia de pelo menos uma raiz, mas pode ser que o intervalo contenha mais de uma raiz como mostra os exemplos na Figura 2.2. Pelo exemplos podemos notar que se f0(x) preserva o sinal em [a, b]e f(a)f(b) < 0, ent˜ao o intervalo cont´´ao podemos afirmar nada sobre em uma unica raiz. Se f(a)f(b) > 0 n˜ a existˆencia ou n˜ao de ra´izes. Exemplo 2.1.2 Da an´alise gr´afica vimos que a fun¸c˜ao f(x)= e¡x - x tem uma raiz em [0, 1] . Tabelando a fun¸c˜ao para valores a partir de zero e espa¸cados de 0.25 temos x 0 0.25 0.5 0.75 1 f(x) 1 0.528 0.106 -0.277 -0.632 Temos que f(0:5)f(0:75) < 0. Logo a raiz pertence ao intervalo [0:5, 0:75]. Note que f0(x)= ¡e¡x ¡1 < 0 8x . IR, isto ´e f. preserva o sinal em [a, b] e com isto podemos concluir que esta raiz ´e ´unica. Devemos observar que o tabelamento ´e uma estrat´egia que completa a an´alise gr´afica. Somente com o tabelamento n˜ao conseguimos determinar se existe outras ra´izes no intervalo ou ainda em que intervalo devemos tabelar a fun¸c˜ao. ˜ CAP´ITULO 2. ZEROS DE FUNC¸OES 13 00.511.522.53-3-2-101234 x1| x2| ab 00.511.522.5300.511.522.5 x | ab f0(x) preserva sinal f0(x) muda de sinal f0(x) muda de sinal f(a)f(b) > 0 Figura 2.2: Exemplos do comportamento de f(x) ˜ CAP´ITULO 2. ZEROS DE FUNC¸OES 14 2.2 Refinamento Nas pr´oximas se¸c˜oes estudaremos os esquemas num´ericos que partindo de uma aproxima¸c˜ao inicial x0, v˜ao gerar uma seq¨uˆencia fxk} que converge para a raiz procurada, isto ´e xk . ! quando k A aproxima¸c˜ao inicial parte do intervalo encontrado na Fase I, Isolamento !1. das Ra´izes, e os termos da seq¨uˆencia s˜ao calculados at´e que a aproxima¸c˜ao tenha atingido uma precis˜ao desejada (crit´erio de parada). 2.3 M´etodo da Bissec¸c˜ao Este m´etodo ´e baseado no Teorema 2.1.1. Seja f(x) uma fun¸c˜ao cont´inua no intervalo [a, b] tal que f(a)f(b) < 0 e seja "> 0 um n´umero dado. A id´eia ´e reduzir a amplitude do intervalo at´e atingir a precis˜ao requerida: b - a<", usando divis˜ao sucessivas do intervalo. Figura 2.3: M´etodo da Bissec¸c˜ao O m´etodo procede da seguinte forma: fa¸ca [a0;b0]=[a, b], x0 = a0 + b0 8> <> f(a0) < 0 f(b0) > 0 8> <> . . (a0;x0) a1 = a0 2 . . : f(x0) > 0 : b1 = x0 ˜ CAP´ITULO 2. ZEROS DE FUNC¸OES 15 8> <> 8> <> f(a1) < 0 . . (x1;b1) a2 = x1 a1 + b1 f(b1) > 0 x1 = . . 2 : : f(x1) < 0 b2 = b1 8> <> 8> <> f(a2) < 0 . . (x2;b2) a3 = x2 a2 + b2 f(b2) > 0 x2 = . . 2 : : f(x2) < 0 b3 = b2 E assim vamos calculando a seq¨uˆencia xk at´e que seja satisfeito o crit´erio de parada bk - ak < ". Este crit´erio garante se tomarmos ¯x . [ak;bk] teremos a garantia que o erro ´e menor que ", isto ´e jx¯- »j= bk - ak Ep, x=(ao+bo)/2; if f(x)*f(ao) > 0, ao=x; else bo=x; end; end; y=(ao+bo)/2; 2.3.1 Estudo da Convergˆencia A convergˆencia ´e bastante intuitiva, como podemos ver na Figura 2.3. Vamos dar uma demonstra¸c˜ao anal´itica atrav´es do seguinte teorema: Teorema 2.3.1 Seja f uma fun¸c˜ao cont´inua em [a, b], onde f(a)f(b) < 0. Ent˜ao o m´etodo da Bissec¸c˜ao gera uma seq¨uˆencia fxk} que converge para a raiz . quando k !1. Prova: O m´etodo gera trˆes seq¨uˆencias: ˜ CAP´ITULO 2. ZEROS DE FUNC¸OES 16 fakg: Seq¨uˆencia n˜ao decrescente e limitada superiormente por b0. Logo tal que lim ak = Ma0 = a1 ·¢¢· a0 )9m . IR k!8 fxkg: Por constru¸c˜ao temos que ak + bk xk = 2 . ak b0 - a0 . 2k . e Como estes valores s˜ao sempre positivos, podemos aplicar a fun¸c˜ao logaritmo, obtendo, log(b0 - a0) - log(") k> log(2) ˜ CAP´ITULO 2. ZEROS DE FUNC¸OES 17 Exemplo 2.3.1 No exemplo 2.1.2 isolamos uma raiz de f(x)= e¡x - x no intervalo [0:5, 0:75]. Usando a precis˜ao e = 10¡8, temos log(0:75 - 0:5) - log(10¡8) k> = 24:575. log(2) Logo ser´a necess´ario no m´inimo 25 itera¸c˜oes para que o m´etodo da Bissec¸c˜ao possa atingir a precis˜ao desejada. 2.4 M´etodo Iterativo Linear (M.I.L.) Seja f(x) cont´inua em [a, b], onde existe uma raiz da equa¸c˜ao f(x) = 0. A estrat´egia deste m´etodo ´e escrever a fun¸c˜ao f de tal forma que f(x)= x - Á(x). Se f(x) = 0, ent˜ao x - Á(x)=0 x = Á(x). Isto ´e, encontrar as ra´izes de f ´e equivalente a achar os pontos fixo da fun¸c˜ao Á. Atrav´es da equa¸c˜ao acima montamos um processo iterativo, onde, dado x0 xn+1 = Á(xn);n =1, 2;::. A fun¸c˜ao f ´e chamada de fun¸c˜ao de itera¸ao e esta n˜ao ´e determinada de forma unica. c˜´ As condi¸c˜oes de convergˆencia s˜ao dadas no teorema abaixo. Teorema 2.4.1 Seja . uma raiz da fun¸c˜ao f isolada no intervalo [a, b]. Seja f uma fun¸c˜ao de itera¸c˜ao da fun¸c˜ao f que satisfaz: 1) f e Á. s˜ao cont´inuas em [a, b], 2) jÁ0(x)j= M< 1 8x . [a, b], 3) x0 . [a, b]. Ent˜ao a seq¨uˆencia fxk} gerada pelo processo iterativo xn+1 = Á(xn) converge para ». Prova: Sendo . uma raiz ent˜ao f(»)=0 . = Á(»), logo . xn+1 = Á(xn) xn+1 - . = Á(xn) - Á(»). . Como f ´e cont´inua e diferenci´avel, pelo Teorema do Valor M´edio temos que existe cn pertencente ao intervalo entre cn e . tal que Á(xn) - Á(»)= Á0(cn)(xn - ») Logo jxn+1 - »| = jÁ0(cn)jjxn - »j= Mjxn - »| ˜ CAP´ITULO 2. ZEROS DE FUNC¸OES 18 Aplicando esta rela¸c˜ao para n - 1;n - 2;, 0 e usando o fato que x0 . [a, b] temos ¢¢· jxn+1 - »j= Mn+1 jx0 - »| Como M< 1, aplicando o limite para n !8 segue que 0 = lim xn+1 - »= lim Mn+1 =0 n!8 j| n!8 jx0 - »| Logo lim xn+1 = . n!8 Observamos que quanto menor for Á0(x)mais r´apida ser´a a convergˆencia. j| Exemplo 2.4.1 Consideremos a fun¸c˜ao f(x)= e¡x¡x, onde existe uma raiz . . [0:5, 0, 75]. Uma forma de escrever f(x)= x - Á(x) ´e considerar Á(x)= e¡x . Verificando as condi¸c˜oes de convergˆencia temos: 1) As fun¸c˜oes Á(x)= e¡x e Á0(x)= ¡e¡x s˜ao cont´inuas em [0:5, 0:75]. 2) A fun¸c˜ao Á. satisfaz max Á0(x)=0:6065::. < 1 (Por que? Ver Nota 1) x2[0:5;0:75] j| 3) Tomando x0 . [0:5, 0:75] teremos garantia de convergˆencia, por exemplo podemos tomar x0 como o ponto m´edio do intervalo 0:5+0:75 x0 = =0:625 2 Assim temos que x1 = Á(x0)= Á(0:625) = 0:53526::. x2 = Á(x1)= Á(0:53526) = 0:58551::. x3 = Á(x2)= Á(0:58551) = 0:55681::. x4 = Á(x3)= Á(0:55681) = 0:57302::. x5 = Á(x4)= Á(0:57302) = 0:56381::. x6 = Á(x5)= Á(0:56381) = 0:56903::. Na Figura 2.4 podemos ver que o comportamento do processo iterativo converge para a raiz. .. . .. . .. . ˜ CAP´ITULO 2. ZEROS DE FUNC¸OES 19 Figura 2.4: M´etodo Iterativo Linear 2.4.1 Crit´erio de Parada Uma quest˜ao ainda est´a em aberto. Qual o xn que fornece uma aproxima¸c˜ao para a raiz, com uma certa precis˜ao e dada. Neste caso podemos usar como crit´erio de parada as seguintes condi¸c˜oes jxn+1 - xnj= e (Erro Absoluto) jxn+1 - xnj= e (Erro Relativo) jxn+1| e vamos tomar xn+1 como aproxima¸c˜ao para a raiz. Se no exemplo anterior tiv´essemos escolhido e =0:006 e o Erro Absoluto ter´iamos jx1 - x0| = j0:53526 - 0:625| =0:08974 >e jx2 - x1| = j0:58551 - 0:53526| =0:05025 >e jx3 - x2| = j0:55681 - 0:58551| =0:02870 >e jx4 - x3| = j0:57302 - 0:55681| =0:01621 >e jx5 - x4| = j0:56381 - 0:57302| =0:00921 >e jx6 - x5| = j0:56903 - 0:56381| =0:00522 e x2 = x1 + e e ¡ ¡ x x 11 ¡ +1 x1 =0:56714 jx2 - x1| =0:0006 Ep, xo=x1; x1=xo-f(xo)/df(xo) end; 2.6 Ordem de Convergˆencia Na se¸c˜ao anterior determinamos o M´etodo de Newton-Raphson que pode ser interpretado como um caso particular do M´etodo Iterativo Linear, onde a convergˆencia ´e mais r´apida. A “ medida” que permite comparar a convergˆencia entre os m´etodos ´e o que chamamos de ordem de convergˆencia, definida por: ˜ CAP´ITULO 2. ZEROS DE FUNC¸OES 23 Defini¸c˜ao 2.6.1 Seja fxn} uma seq¨uˆencia que converge para um n´umero . e seja ek = xk ¡. o erro na itera¸c˜ao k. Se lim jek+1 p| = C p> 1 C> 0 k!8 jekj dizemos que a seq¨uˆencia converge com ordem p e com constante assint´otica C. Como a seq¨uˆencia converge, para valores de k suficientemente grande temos jek+1j˜ Cjekjp, com jek| < 1 Assim quanto maior for o valor de p, menor ser´a o erro ek+1. Quando p = 1 dizemos que o jj m´etodo tˆem convergˆencia linear. Se p = 2 dizemos que a convergˆencia ´e quadr´atica. Primeiramente vamos determinar a ordem de convergˆencia do M.I.L. Sendo a seq¨uˆencia fxn} gerada por xk+1 = Á(xk);k =0, 1, 2;::. e que . = Á(») temos xk+1 - . = Á(xk) - Á(»)= Á0(ck)(xk - »), onde a ´ultima igualdade ´e conseq¨uˆencia do Teorema do Valor M´edio e ck ´e um n´umero entre xk e ». Logo segue xk+1 - . = Á0(ck) ek+1 = Á0(ck) xk - . . ek Aplicando o m´odulo e calculando o limite quando k tende ao infinito temos lim jek+1| = lim Á0(ck)= Á0(»)= C ek k!8 j| k!8 jjj| Portanto temos que o M.I.L. tˆem ordem de convergˆencia p = 1. No caso do M´etodo de Newton-Raphson temos que a seq¨uˆencia ´e gerada pelo processo iterativo f(xn) xn+1 = xn - f(xn) Subtraindo . de cada lado temos f(xn) f(xn) xn+1 - . = xn - . - f(xn en+1 = en - f(xn) (2.2) 0) . 0 Atrav´es da f´ormula de Taylor da fun¸c˜ao f no ponto xn temos f(x)= f(xn)+ f0(xn)(x - xn)+ f00(cn) (x - xn)2 cn . [x, xn] 2 Que calculada em x = . fornece 0= f(»)= f(xn)+ f0(xn)(. - xn)+ f00(cn) (. - xn)2 2 CAP´ITULO 2. ZEROS DE FUNC¸OES˜24 Dividindo por f0(xn) e fazendo en = xn - . segue que f(xn) f00(cn) 2 f0(xn) = en - 2f0(xn) en Substituindo em (2.2) obtemos en+1 f00(cn) = en 2 2f0(xn) Finalmente aplicamos o m´odulo e calculamos o limite quando k tende ao infinito obtendo lim jen+1| = lim . f00(cn) . = . f00(») . = 1 Á00(»)= C 2 k!8 jenjk!8 2f0(xn)2f0(»)2j| Portanto temos que o M´etodo de Newton-Raphson tˆem ordem de convergˆencia p = 2. 2.7 Observa¸c˜oes Finais Neste cap´itulo vimos trˆes m´etodos diferentes para resolver equa¸c˜oes da forma f(x) = 0. Faremos um breve coment´ario das vantagens e desvantagens de cada m´etodo. No M´etodo da bissec¸c˜ao vimos que o n´umero de itera¸c˜oes depende apenas do intervalo inicial [a0;b0] Logo este pode ser aplicado a qualquer fun¸c˜ao f(x) que satisfaz f(a)f(b) < 0. N˜ao importa o quanto f(x) seja complicada. A desvantagem ´e que tem uma convergˆencia lenta. Na pr´atica ele ´e usado para refinar o intervalo que cont´em a raiz. Aplicamos o m´etodo em um n´umero fixo de itera¸c˜oes. Em geral o M.I.L. ´e mais r´apido que o M´etodo da Bissec¸c˜ao. Usa menos opera¸c˜oes por cada itera¸c˜ao. Pode encontrar ra´izes em intervalos onde f(a)f(b) > 0 . A dificuldade ´e encontrar a fun¸c˜ao de itera¸c˜ao f que seja convergente. O M´etodo de Newton-Raphson tˆem convergˆencia quadr´atica. Por´em este necessita da avalia¸c˜ao da fun¸c˜ao e sua derivada em cada ponto xn. Pode ocorrer de termos uma raiz isolada num intervalo [a, b] e o m´etodo acabe convergindo para uma outra raiz que n˜ao ¯ pertence a [a, b]. Isto ocorre porque temos que tomar x0 . [¯a, b] . [a, b]. Na pr´atica tomamos x0 como ponto m´edio do intervalo, isto ´e a + b x0 = 2 Nota 1 Em muitas situa¸c˜oes vamos necessitar de calcular o m´aximo do m´odulo de uma fun¸c˜ao restrita a um intervalo, isto ´e max f(x). x2[a;b] jj Uma forma pr´atica para este c´alculo ´e seguir os passos: 1: Calcula-se os valores da fun¸c˜ao nos extremos do intervalo, f(a)e f(b). jjjj ˜ CAP´ITULO 2. ZEROS DE FUNC¸OES 25 2: Verifique se a fun¸c˜ao n˜ao possui ponto critico no intervalo, ou seja, achamos os valores de xk tal que f0(xk)=0 e xk . [a, b] 3: Tomamos como o valor m´aximo o maxfjf(a)j, jf(b)j, jf(xk)j} 2.8 Exerc´icios Exerc´icio 2.1 Localize graficamente e dˆe intervalos de amplitude 0.5 que contenha as ra´izes das equa¸c˜oes a) ln(x)+2x =0 b) ex sen(x)=0 c) ln(x) - 2x = ¡2 x 2 e- x d) 2 cos(x) - =0 e) 3 ln(x) - f) (5 - x)ex =1 22 Exerc´icio 2.2 Utilize o M´etodo da Bissec¸c˜ao e aproxime a menor raiz em m´odulo com erro relativo menor que 10¡1 para as equa¸c˜oes a) e b) do exerc´icio anterior. Exerc´icio 2.3 Utilize o M´etodo Iterativo Linear e aproxime a menor raiz em m´odulo com erro relativo menor que 10¡2 para as equa¸c˜oes c) e d) do exerc´icio anterior. Exerc´icio 2.4 Utilize o M´etodo de Newton-Rapshon e aproxime a menor raiz em m´odulo com erro relativo menor que 10¡3 para as equa¸c˜oes d) e f) do exerc´icio anterior. Exerc´icio 2.5 Achar a raiz p-´esima de um n´umero positivo a ´e equivalente a achar a raiz p positiva da equa¸c˜ao pa = x. a) Encontre um intervalo que depende do valor de a e que contenha a raiz. b) Verifique se a fun¸c˜ao de itera¸c˜ao Á(x)= a=xp¡1 satisfaz os crit´erios de convergˆencia do M´etodo Iterativo Linear. c) Verifique que o processo iterativo gerado pelo M.N.R. ´e dado por 1 ·a . xn+1 = p (p - 1)xn + xnp¡1 d) Implemente um programa em Matlab que execute o processo iterativo dado em c). Exerc´icio 2.6 Dada a fun¸c˜ao f(x)= ex - 4x2 . a) Isole as ra´izes da fun¸c˜ao f(x). b) Verifique que as fun¸c˜oes abaixo s˜ao fun¸c˜ao de itera¸c˜ao de f e verifique se satisfazem o crit´erio de convergˆencia do M.I.L. para a raiz positiva. Á1(x)= 1 ex=2 Á2(x) = ln(4x 2) 2 ˜ CAP´ITULO 2. ZEROS DE FUNC¸OES 26 c) Tomando x0 =0:6 e e =0:01, aplique o M.I.L. para encontrar uma aproxima¸c˜ao para a raiz positiva, usando uma fun¸c˜ao de itera¸c˜ao que satisfa¸ca os crit´erios de convergˆencia Exerc´icio 2.7 Prove o Teorema 2.5.1. Exerc´icio 2.8 A fun¸c˜ao f(x)= sen(cos(p3x)) tem uma raiz no intervalo [0:7, 0:9]. Encontre uma aproxima¸c˜ao com e =0:07, escolhendo entre os m´etodos num´ericos estudados o mais adequado. Justifique sua resposta. Cap´itulo 3 Sistemas Lineares A resolu¸c˜ao de sistemas lineares surge em diversas ´areas do conhecimento. O caso geral em que o sistema linear envolve m equa¸c˜oes com n inc´ognitas, o sistema pode apresentar uma unica solu¸´c˜ao, infinitas solu¸c˜oes ou n˜ao admitir solu¸c˜ao. Este tipo de problema ´e tratado na ´ Algebra Linear usando o processo de escalonamento. Neste cap´itulo vamos analisar esquemas num´ericos para solu¸c˜oes de sistemas lineares de n equa¸c˜oes com n inc´ognitas , isto ´e 8 . a1;1x1 + a1;2x2 + a1;3x3 a1;nxn = b1 ¢¢· . a2;1x1 + a2;2x2 + a2;3x3 ¢¢· a2;nxn = b2 a3;1x1 + a3;2x2 + a3;3x3 a3;nxn = b3 ¢¢· ... .. . ... .. ... .. . an;1x1 + an;2x2 + an;3x3 an;nxn = bn ¢¢· onde aij s˜ao os coeficientes, xj s˜ao as inc´ognitas e os bj s˜ao as constantes. Este sistema pode ser escrito na forma matricial Ax = b com A . IRn£n e x, b . IRn . Analisaremos duas classes de esquemas num´ericos: M´etodos Diretos e M´etodos Iterativos. 3.1 M´etodos Diretos Os M´etodos Diretos s˜ao aqueles que ap´os um n´umero finito de opera¸c˜oes fornecem a solu¸c˜ao exata do sistema, a menos dos erros de arredondamentos. Estes m´etodos s˜ao baseados no ´ processo de escalonamento estudado em Algebra Linear. S˜ao eficientes para sistemas de pequeno porte (n˜ao mais que 50 equa¸c˜oes ) e para sistemas de bandas, como por exemplo sistemas tridiagonais ( ver Ex. 3.3 ). Primeiramente vamos considerar os sistemas lineares triangulares. 27 CAP´ITULO 3. SISTEMAS LINEARES 28 3.1.1 Sistema Triangular Superior Um Sistema Triangular Superior ´e aquele em que a matriz associada ao sistema ´e uma matriz triangular superior, isto ´e ai;j = 0 para i>j. 8. <. . a1;1x1 + a1;2x2 + a1;3x3 a1;nxn = b1 ¢¢· a2;2x2 + a2;3x3 a2;nxn = b2 ¢¢· a3;3x3 a3;nxn = b3 ¢¢· .. .. .. an;nxn = bn Este sistema admite uma ´c˜= 0 para i =1, 2;:::;n, sendo, unica solu¸ao se aii bn xn = an;n 1 xn¡1 = an¡1;n¡1 (bn¡1 - an¡1;nxn) 1 . xn¡2 = an¡2;n¡2 (bn¡2 - an¡2;n¡1xn¡1 - an¡2;nxn) .. .. .. n 0@ 1. 1 xk = bk - ak;jxj ak;k j=k+1 . . . . . . e assim sucessivamente. Com isto obtemos o esquema num´erico para solu¸c˜ao de sistema triangular superior dado pelo algoritmo abaixo Algoritmo: Retro-Solu¸c˜ao Input: Matriz triangular superior A . IRn£n e b . IRn xn bn=an;n . 1. 0@ Para k = n - 1;n - 2;::. 1, fa¸ca: n bk - 1 xk . ak;k ak;jxj j=k+1 fim para Output: x . IRn : solu¸c˜ao do sistema CAP´ITULO 3. SISTEMAS LINEARES 29 3.1.2 M´etodo de Elimina¸c˜ao de Gauss Dois sistemas lineares s˜ao ditos ser equivalentes se estes tem a mesma solu¸c˜ao. A estrat´egia do M´etodo de Elimina¸c˜ao de Gauss ´e transformar um sistema linear Ax = b em um sistema triangular superior equivalente Sx = b˜ cuja a solu¸c˜ao ´e facilmente obtida pela Retro-Solu¸c˜ao. Esta transforma¸c˜ao ´e realizada atrav´es das opera¸c˜oes elementares (I) Trocar duas equa¸c˜oes. (II) Multiplicar uma equa¸c˜ao por uma constante n˜ao nula. (III) Adicionar a uma equa¸c˜ao uma outra multiplicada por uma constante n˜ao nula. Aplicando qualquer seq¨uˆencia dessas opera¸c˜oes elementares num sistema Ax = b obtemos um novo sistema ˜b˜ de tal forma que estes ser˜etodo Ax = ao equivalentes. Para descrever o M´ de Elimina¸c˜ao de Gauss vamos considerar o sistema linear 8. <. . a1;1x1 + a1;2x2 + a1;3x3 · · · a1;nxn = b1 a2;1x1 + a2;2x2 + a2;3x3 · · · a2;nxn = b2 a3;1x1 + a3;2x2 + a3;3x3 · · · a3;nxn = b3 , . . . . . . . . . . . . . . . an;1x1 + an;2x2 + an;3x3 · · · an;nxn = bn onde det(A) = e, o sistema admite uma unica solu¸c˜ao. 60, isto ´´Um sistema linear pode ser representado na forma de matriz estendida (A0b0), ou seja 0. . | (0) (0) (0) (0)a1;1 a1;2 a1;3 a1;n (0) (0) (0) ¢¢· (0)a2;1 a2;2 a2;3 a2;n (0) (0) (0) ¢¢· (0)a3;1 a3;2 a3;3 a3;n ¢¢· ... . ... . ... . (0) (0) (0) (0)an;1 an;2 an;3 an;n ¢¢· (0)b1 (0)b2 (0)b3 . . . b(0) n 1. . onde o ´indice superior indica a etapa do processo. (0) Etapa 1: Eliminar a inc´ognita x1 das equa¸c˜oes k =2, 3;:::;n. Sendo a1;1 = 0, usaremos a opera¸c˜ao elementar (III) e subtra´imos da linha k a primeira linha multiplicada por (0)ak;1 mk;1 = (0) . a1;1 (0) Os elementos mk;1 s˜ao chamados de multiplicadores e o elemento a1;1 ´e chamado de (0) pivˆo da Etapa 1. Indicando a linha k da matriz entendida por Lk esta etapa se resume em L(1) 1 = L(0) 1 L(1) k = L(0) k - mk;1L(0) 1 , k = 2, 3, . . . , n CAP´ITULO 3. SISTEMAS LINEARES 30 Ao final desta etapa teremos a matriz 1. . (1)b1 (1)b2 (1)b3 . . . b(1) n (1) (1) (1) (1)a1;1 a1;2 a1;3 a1;n (1) (1) ¢¢· (1) 0 a2;2 a2;3 a2;n (1) (1) ¢¢· (1) 0 a3;2 a3;3 a3;n ¢¢· ... . ... . ... . (1) (1) (1)0 an;2 an;3 a ¢¢· n;n 0. . que representa um sistema linear equivalente ao sistema original, onde a inc´ognita x1 foi eliminada das equa¸c˜oes k =2, 3;:::;n. (1) Etapa 2: Eliminar a inc´ognita x2 das equa¸c˜oes k =3, 4;:::;n. Supondo que a2;2 = 0, vamos tomar este elemento como pivˆo desta etapa e desta forma os multiplicadores s˜ao dados por (1)ak;2 mk;2 = (1)a2;2 A elimina¸c˜ao se procede da seguinte forma: L(2) 1 = L(1) 1 L(2) 2 = L(1) 2 L(2) k = L(1) k - mk;2L(1) 2 , k = 3, 4, . . . , n obtendo ao final da etapa a matriz Com procedimentos an´alogos ao das etapas 1 e 2 podemos eliminar as inc´ognitas xk das equa¸c˜oes k +1;k +2;:::;n e ao final de n - 1 etapas teremos a matriz Esta matriz representa um sistema triangular superior equivalente ao sistema original. Logo a solu¸c˜ao deste sistema, obtido pela Retro-Solu¸c˜ao, ´e solu¸c˜ao do sistema original. 1. . b1(n¡1) b2(n¡1) b3(n¡1) . . . b(n¡1) n a1(n ;1¡1) a1(n ;2¡1) a1(n ;3¡1) a1(n ;n¡1) (n¡1) (n¡1) ¢¢· (n¡1)0 a2;2 a2;3 a2;n (n¡1) ¢¢· (n¡1)00 a3;3 a3;n ¢¢· ... . ... . ... . 000 a(n¡1) ¢¢· n;n 1. . (2)b1 (2)b2 (2)b3 . . . b(2) n a(2) 1;1 a(2) 1;2 a(2) 1;3 · · · a(2) 1;n 0 a(2) 2;2 a(2) 2;3 · · · a(2) 2;n 0 0 a(2) 3;3 · · · a(2) 3;n . . . . . . . . . . . . 0 0 a(2) n;3 · · · a(2) n;n 0. . 0. . CAP´ITULO 3. SISTEMAS LINEARES 31 Algoritmo: M´etodo de Elimina¸c˜ao de Gauss Input: Matriz A e vetor b . IRn Elimina¸c˜ao: Para k =1, 2;:::;n - 1, fa¸ca: Para i = k +1;:::;n, fa¸ca: aij m . ak;k Para j = k +1;:::;n, fa¸ca: aij aij - m * akj à fim para bi bi - m * bk . fim para fim para Retro-Solu¸c˜ao: xn bn=an;n . Para k = n - 1;n - 2;::. 1, fa¸ca: 1n A bk - 0@ 1 xk ak;jxj . ak;k j=k+1 fim para Output: x . IRn : solu¸c˜ao do sistema Exemplo 3.1.1 Vamos considerar o sistema linear abaixo 8>: <> 3x1 +2x2 - x3 =1 7x1 - x2 - x3 = ¡2 x1 + x3 =1 Escrevendo na forma de matriz estendida teremos . B. 32 ¡1 1 7 ¡1 ¡1 ¡2 101 1 . C. Etapa 1: Eliminar x1 das linhas 2 e 3. (1) (0) L1 = L1 (0) (1) (0) (0) a21 7 L2 = L2 - m2;1L1 , onde m2;1 == (0)a1;1 3 (0) (1) (0) (0) a31 1 L3 = L3 - m3;1L1 , onde m3;1 == (0)a1;1 3 CAP´ITULO 3. SISTEMAS LINEARES 32 e com isto obtemos a matriz . B. 32 ¡1 1 0 ¡17=34=3 ¡13=3 0 ¡2=34=3 12=3 . C. Etapa 2: Eliminar x2 da linha 3. L(2) 1 = L(1) 1 L(2) 2 = L(1) 2 L(2) 3 = L(1) 3 - m3;2L(1) 2 , onde m3;2 = a(01) 32 a(1) 2;2 = 2 17 obtendo assim a matriz . B. 32 ¡1 1 0 ¡17=34=3 ¡13=3 0 020=17 20=17 . C. Retro-Solu¸c˜ao: Encontrar a solu¸c˜ao do sistema triangular superior. b3 x3 = =1 a3;3 1 x2 =(b2 - a2;3x3)=1 a2;2 1 x1 =(b1 - a1;2x2 - a1;3x3)=0 a1;1 Logo a solu¸c˜ao do sistema ´e dada por x = (0, 1, 1)T . A solu¸c˜ao encontrada ´e a solu¸c˜ao exata, pois mantivemos os n´umeros resultantes na forma de fra¸c˜ao. Por´em m´aquinas digitais representam estes n´umeros na forma de ponto flutuante finita e erros de arredondamento podem ocorrer. Em sistemas lineares de grande porte estes erros v˜ao se acumulando e prejudicando a solu¸c˜ao do sistema. 3.1.3 Pivotamento Parcial Em cada etapa k da elimina¸c˜ao temos o c´alculo do multiplicador (k¡1)ak;j mk;j = (k¡1) . ak;k Se o pivˆo ja(k¡1)j. 1, ou seja este ´e pr´oximo de zero teremos problemas com os erros de k;k arredondamento, pois operar n´umeros de grandezas muito diferentes aumenta os erros ( ver Ex. 3.4). A estrat´egia de pivotamento parcial ´e baseada na opera¸c˜ao elementar (I). No in´icio de cada etapa k escolhemos como pivˆo o elemento de maior m´odulo entre os coeficientes aik k¡1 para i = k, k +1;:::;n. CAP´ITULO 3. SISTEMAS LINEARES 33 Exemplo 3.1.2 Vamos considerar o sistema linear, representado pela matriz extendida . B. 1 3 2 3 2 1 2 4 ¡3 2 1 ¡2 . C. Etapa 1: Escolha do pivˆo max ai;1= a3;11·i·3 jjj| Trocamos a linha L1 com a linha L3, obtendo, . B. ¡321 ¡2 212 4 132 3 . C. Eliminar x1 das linhas 2 e 3. L(1) 1 L(1) 2 L(1) 3 = = = L(0) 1 L(0) 2 - m2;1L(0) 1 , L(0) 3 - m3;1L(0) 1 , onde onde m2;1 = 2 3 m3;1 = 1 3 e com isto obtemos a matriz 3- @0 . 21 ¡2 07=3 8=3 8=3 . CA 0 11=37=3 7=3 Etapa 2: Escolha do pivˆo max ai;1= a3;22·i·3 jjj| Trocamos a linha L2 com a linha L3, obtendo, . B. ¡3 21 ¡2 0 11=37=3 7=3 07=38=3 8=3 . C. Eliminar x2 da linha 3. L(2) 1 L(2) 2 L(2) 3 = = = L(1) 1 L(1) 2 L(1) 3 - m3;2L(1) 2 , onde m3;2 = ¡7 11 obtendo assim a matriz . B. ¡32 1 ¡2 0 11=38=3 8=3 0 013=11 13=11 . C. CAP´ITULO 3. SISTEMAS LINEARES 34 Retro-Solu¸c˜ao: Encontrar a solu¸c˜ao do sistema triangular superior. b3 x3 = =1 a3;3 1 x2 =(b2 - a2;3x3)=0 a2;2 1 1C x1 =(b1 - a1;2x2 - a1;3x3)=1 a1;1 Logo a solu¸c˜ao do sistema ´e dada por x = (1, 0, 1)T . Na pr´atica a troca de linhas n˜ao ´e realizada. O controle ´e feito por um vetor de inteiros n-dimensional, onde inicialmente na posi¸c˜ao k est´a armazenado k, ou seja trc = [1, 2, . . . , s, . . . , n]. Se, por exemplo, trocamos a linha 1 pela linha s o vetor passa a ser trc =[s, 2;:::, 1;:::;n]. 3.1.4 C´alculo da Matriz Inversa Vamos supor que desejamos resolver os sistemas lineares Ax = b1 , Ax = b2 ;::. Ax = bk , onde a matriz A ´e a mesma para todos os sistemas. A matriz triangular superior, resultante do processo de elimina¸c˜ao, n˜ao depende do vetor b e portanto ser´a a mesma em qualquer um dos sistemas. Assim podemos resolver estes sistemas num ´c˜ unico processo de elimina¸ao usando a matriz estendida (Ab1b2::. bk) e aplicando a Retro-Solu¸c˜ao para cada vetor bk . jj| j O C´alculo da inversa de uma matriz ´e um caso particular do esquema acima. A inversa de uma matriz A . Rn£n, denotada por A¡1, ´e uma matriz n × n tal que AA¡1 = A¡1 A = I Como exemplo vamos considerar uma matriz A de dimens˜ao 3 × 3 . , . B. 1C . 1C . , 1C . = a1;1 a1;2 a1;3 a2;1 a2;2 a2;3 a3;1 a3;2 a3;3 cuja a inversa A¡1 ´e dada por . B. x1;1 x1;2 x1;3 x2;1 x2;2 x2;3 x3;1 x3;2 x3;3 logo temos que . B. . B. . C. . B. 100 a1;1 a1;2 a1;3 x1;1 x1;2 x1;3 010 a2;1 a2;2 a2;3 x2;1 x2;2 x2;3 a3;1 a3;2 a3;3 x3;1 x3;2 x3;3 001 CAP´ITULO 3. SISTEMAS LINEARES 35 Portanto cada coluna k da inversa da matriz A ´e solu¸c˜ao de um sistema linear, onde o vetor dos termos independentes ´e a k-´esima coluna da matriz identidade, isto ´e 1C =A. , 1C =A. , 1C =A. , 1C 1C1C . B. a1;1 a1;2 a1;3 a2;1 a2;2 a2;3 . B. . C. x1;1 x2;1 . B. 1 0 0 a3;1 a3;2 a3;3 x3;1 . B. a1;1 a1;2 a1;3 a2;1 a2;2 a2;3 . B. . C. x1;2 x2;2 . B. 0 1 0 a3;1 a3;2 a3;3 x3;2 . B. a1;1 a1;2 a1;3 a2;1 a2;2 a2;3 . B. . C. x1;3 x2;3 . B. 0 0 1 a3;1 a3;2 a3;3 x3;3 Em resumo, se temos uma matriz n × n, podemos achar a inversa resolvendo n sistemas lineares, representados pela matriz estendida (Ab1b2::. bn), onde os vetores bk s˜ao os jj| j 1C vetores unit´arios ( 1 na posi¸c˜ao k e zeros nas demais posi¸c˜oes). Exemplo 3.1.3 Vamos achar a inversa da matriz abaixo, usando o m´etodo de Elina¸c˜ao de . , Gauss. . B. 41 ¡6 32 ¡6 ¡5 31 1C Para o processo de elimina¸c˜ao consideremos a matriz estendida . , . B. 41 ¡6 1 0 0 32 ¡6 0 1 0 ¡5 31 0 0 1 1C Etapa 1: Eliminar x1 das linhas 2 e 3. (1) (0) L1 = L1 (1) (0) (0) 3 L2 = L2 - m2;1L1 , onde m2;1 = 4 (1) (0) (0) 3 L3 = L3 - m3;1L1 , onde m3;1 = 4 e com isto obtemos a matriz . , . B. 41 ¡6 1 0 0 05=4 ¡3=2 1 0 ¡3=4 ¡1=2 ¡3=4 01=4 0 1 CAP´ITULO 3. SISTEMAS LINEARES 36 Etapa 2: Eliminar x2 da linha 3. (2) (1) L1 = L1 (2) (1) L2 = L2 (2) (1) (1) 1 L3 = L3 - m3;2L2 , onde m3;2 = 5 obtendo assim a matriz . C. , . B. 41 ¡6 1 0 0 05=4 ¡3=2 1 0 ¡3=4 ¡1=5 ¡3=5 00 1 ¡1=5 Retro-Solu¸c˜ao: Encontrar a solu¸c˜ao dos sistemas triangulares superior. Primeira coluna da inversa x3;1 = b1 3 a3;3 = 3 x2;1 = 1 a2;2 (b1 2 - a2;3x3) = 3 x1;1 = 1 a1;1 (b1 1 - a1;2x2 - a1;3x3) = 4 Segunda coluna da inversa x3;2 = b2 3 a3;3 = 1 x2;2 = 1 a2;2 (b2 2 - a2;3x3) = 2 x1;2 = 1 a1;1 (b2 1 - a1;2x2 - a1;3x3) = 1 Terceira coluna da inversa x3;3 = b3 3 a3;3 = ¡5 x2;3 = 1 a2;2 (b3 2 - a2;3x3) = ¡6 x1;3 = 1 a1;1 (b3 1 - a1;2x2 - a1;3x3) = ¡6 1C Logo a matriz inversa s ´e dada por . , . B. 41 ¡6 32 ¡6 ¡5 31 CAP´ITULO 3. SISTEMAS LINEARES 37 No exemplo acima temos que a inversa da matriz A ´e a pr´opria A. Este tipo de matriz ´e usado como matriz teste para verificar a eficiˆencia dos m´etodos num´ericos. Abaixo apresentamos uma implementa¸c˜ao do M´etodo de Elimina¸c˜ao de Gauss em MatLab que resolve k sistemas, onde a matriz A ´e comum a todos. % Disciplina de C\'{a}lculo Num\'{e}rico -Prof. J. E. Castilho % M\'{e}todo de elimina\c{c}\~{a}o de Gauss % Este procedimento resolve k-sistemas lineares, onde a % matriz A de dimens\~{a}o n x n eh comum a todos. % Os par\^{a}metros devem ser passados da forma % x=EGauss(A,b1,b2,b3,...,bk) % o resultado eh uma matriz x de dimens\~{a}o n x k onde a % coluna j armazena a solu\c{c}\~{a}o do sistema Ax=bj % Dados A: matriz do sistema % varargin: lista dos vetores dos termos independentes function x=EGauss(A,varargin) b=[varargin{:}]; db=size(b); % Esta parte verifica se o sistema eh quadrado da=size(A); n=da(1); if n ~=da(2), disp('??? A matriz deve ser quadrada'); break; end; % Esta parte verifica se a dimens\~{a}o do vetor b % esta de acordo com a dimens\~{a}o do sistema if n ~=db(1), disp('??? Erro na dimens\~{a}o de b'); break; end; % Cria matriz estendida Ax=[A,b]; % Fase da elimina\c{c}\~{a}o for k=1:(n-1) for i=k+1:n if abs(Ax(k,k)) < 10^(-16), disp('??? Piv\^{o} Numericamente Nulo'); break; end; m=Ax(i,k)/Ax(k,k); for j=k+1:da(2) + db(2) Ax(i,j) = Ax(i,j)-m*Ax(k,j); end; CAP´ITULO 3. SISTEMAS LINEARES 38 end; end; % Fase da Retro solu\c{c}\~{a}o if abs(Ax(n,n)) < 10^(-16), disp('??? Piv\^{o} Numericamente Nulo'); break; end; for m=1:db(2) x(n,m) = Ax(n,n+m)/Ax(n,n); end; for k=(n-1):-1:1 for m=1:db(2) som=0; for j=k+1:n som=som+Ax(k,j)*x(j,m); end; x(k,m) = (Ax(k,n+m)-som)/Ax(k,k); end; end; 3.2 M´etodos Iterativos Vamos considerar um sistema linear Ax = b, onde A . IRn£n e x, b . IRn . Os m´etodos iterativos seguem um esquema semelhante aos m´etodos para o c´alculo de zeros de fun¸c˜oes. O sistema linear ´e escrito na forma x = Cx + g, onde g . IRn e C . IRn£n ´e chamada de matriz de itera¸c˜ao. Assim montamos o processo iterativo: Dado x0 x k+1 = Cxk + g. Sendo um processo iterativo, necessitamos de um crit´erio de parada. E para isto temos que ter uma medida entre as aproxima¸c˜oes xk+1 e xk . Para isto vamos usar o conceito de norma de matrizes, definida abaixo Defini¸c˜ao 3.2.1 Uma norma em IRn£m ´e uma aplica¸c˜£m IR que satisfaz as ao j| · j| : IRn! seguintes propriedades: (M1) jjAj| = 0 e jjAj| =0 . A =0, 8A . IRn£m . CAP´ITULO 3. SISTEMAS LINEARES 39 (M2) jj®Aj| = j®| jjAjj, 8a . IR e 8A . IRn£m . (M3) jjA + Bj| = jjAj| + jjBjj, 8A, B . IRn£m . As normas matriciais mais usadas s˜ao . n 1·j·m m . i=1 . jaij| Norma do M´aximo das Coluna jjAjj1 = max 8. . 9. ; jjAjj8 Norma do M´aximo das Linha = max jaijj 1·i·n j=1 0@ 2 1. 1=2 nm . j=1 . i=1 jaijj jjAjj2 = Norma Euclidiana A norma vetorial pode ser vista como um caso particular da norma matricial, onde um vetor x . IRn ´e equivalente a uma matriz de ordem n × 1. Com isto temos as normas de vetores dadas por jjxjj1 = n. i=1 jxi| Norma da Soma . n. = max xiNorma do M´aximojjxjj8 1·i·n j| !1=2 2 Norma Euclidiana jjxjj2 = jxij i=1 O conceito de norma nos permite definir convergˆencia de uma seq¨uˆencia de vetores fxkg. Dizemos que xk x¯se . lim k x¯j| =0 k!8 jjx . Com isto podemos definir os crit´erios de parada: Dado um "> 0 jjx k+1 - x k j| = e Erro Absoluto k+1k jjxk- xj| = e Erro Relativo jjxj| Teste do Res´iduo jjb - Axk j| = e Al´em disso, as normas j| · jj1, j| · jj2 e j| · jj8 satisfazem as seguintes propriedades: (M4) jjAxj| ·jjAj| jjxj| (M5) jjABj| = jjAj| jjBj| CAP´ITULO 3. SISTEMAS LINEARES 40 3.2.1 Crit´erio de Convergˆencia Dependendo da forma da matriz C a seq¨uˆencia gerada pelo processo iterativo pode ou n˜ao convergir para a solu¸c˜ao do sistema. Seja x¯a solu¸c˜ao do sistema Ax = b, logo x¯satisfaz x¯= C¯x + g. Com isto temos que x k+1 - x¯= C(x k+1 - x¯) Sendo o erro em cada itera¸c˜ao dado por ek = xk - x¯e usando as propriedades de norma (M4) e (M5) segue que jje k j| = jjCj| jj2 e kk ¡ ¡ 12j | = jjCj| jje jj .. .. .. = jjCj| k jje 0 j| Logo a seq¨uˆencia fxk} converge para a solu¸c˜ao do sistema x¯se lim k lim k 0 j| =0 k!8 jje j| = k!8 jjCj| jje e isto ocorre se e somente se a matriz C satisfaz a condi¸c˜ao jjCj| < 1 Note que o crit´erio de convergˆencia n˜ao depende do vetor inicial x0 . A escolha de x0 influˆencia no n´umero de itera¸c˜oes necess´arias para atingir a precis˜ao desejada. Quanto menor for jjx0 - x¯j| menos itera¸c˜oes ser˜ao necess´arias. 3.2.2 M´etodo Iterativo de Gauss-Jacobi Vamos considerar o sistema linear Ax = b dado por 8a1;1x1 + a1;2x2 + a1;3x3 a1;nxn = b1 . ¢¢· a2;1x1 + a2;2x2 + a2;3x3 a2;nxn = b2. a3;1x1 + a3;2x2 + a3;3x3 ¢¢· a3;nxn = b3 ¢¢· , ... .. . .. .. .. .. .. . an;1x1 + an;2x2 + an;3x3 an;nxn = bn ¢¢· CAP´ITULO 3. SISTEMAS LINEARES 41 onde os aii = 0 para i =1, 2;:::;n. Em cada equa¸c˜ao i podemos isolar a inc´ognita xi obtendo as seguintes rela¸c˜oes 1 x1 = a1;1 (b1 - a1;2x2 - a1;3x3 - · · · - a1;nxn) 1 x2 = a2;2 (b2 - a2;1x1 - a2;3x3 - · · · - a2;nxn) 1 x3 = a3;3 (b3 - a3;1x1 - a3;2x2 - · · · - a3;nxn) . . . . . . . . . 1 xn = an;n (bn - an;1x1 - an;2x2 - · · · - an;n¡1xn¡1) Na forma matricial estas equa¸c˜oes s˜ao equivalentes `a 0. 1. x1 1. 0. b1 a1;1 1. 0. x1 0. 1. a1;2 a1;3 a1;n 0 - - ¢¢· - a1;1 a1;1 a1;1 x2 a2;1 a2;3 a2;n x2 b2 ¡a2;2 0 ¡a2;2 ¢¢· - a2;2 a2;2 x3 x3 = + (3.1) .. . a3;1 a3;2 . .0 a3;n . b3 . a3;3 a3;3 a3;3 . . - - ¢¢· - . a3;3 .. . . ... ..... . . ..... . . . .. an;1 an;2 an;3 @ @A 0 @A @A bn an;n . . . - - - ¢¢· xn xn an;n an;n an;n Desta forma temos o sistema linear na forma x = Cx + g e assim montamos o processo iterativo conhecido como M´etodo Iterativo de Gauss Jacobi: Dado x0 (k+1) (k)(k)(k)x1 = a 1 1;1 ³b1 - a1;2x2 - a1;3x3 ¡¢¢¢- a1;nxn (k) n x2(k+1) = a 1 2;2 ³b2 - a2;1x1(k) - a2;3x3(k) ¡¢¢¢- a2;nx x3(k+1) = 1 ³b3 - a3;1x1(k) - a3;2x2(k) n (k) a3;3 ¡¢¢¢- a3;nx . .. . .. . .. (3.2) (k+1) 1 (k)(k)(k)xn = an;n ³bn - an;1x1 - an;2x2 ¡¢¢¢- an;n¡1xn¡1 CAP´ITULO 3. SISTEMAS LINEARES 42 Algoritmo: M´etodo Iterativo de Gauss-Jacobi Input: Matriz A . IRn£n b, x0 . IRn e "> 0 Enquanto jjxk+1 - xkj| >e fa¸ca: Para s =1, 2;:::;n, fa¸ca: n 0@ 1. 1 (k) (k+1) bs ¡ x as;jx j s . as;s j=1;j=s fim para fim enquanto Output: x . IRn : solu¸c˜ao do sistema 3.2.3 Crit´erio das Linhas Como crit´erio de convergˆencia, vimos que a matriz de itera¸c˜ao C deve satisfazer a condi¸c˜ao jjCj| < 1. Usando a Norma do M´aximo das Linhas sobre a matriz C em (3.2) temos o seguinte crit´erio de convergˆencia para o M´etodo de Gauss-Jacobi Teorema 3.2.1 (Crit´erio das Linhas) Dado o sistema linear Ax = b. Seja os ®k de tal forma que: n 1 ®k jak;j| < 1 para k =1, 2;:::;n = j=1;j=k jak;k| 6 Ent˜ao o M´etodo de Gauss-Jacobi gera uma seq¨uˆencia fxk} que converge para a solu¸c˜ao do sistema. Este crit´erio fornece uma condi¸c˜ao necess´aria, mas n˜ao suficiente. Isto ´e, o crit´erio pode n˜ao ser satisfeito e o m´etodo ainda pode convergir. Outros crit´erios podem ser obtidos usando outras normas. Por exemplo, podemos obter o chamado crit´erio das colunas aplicando a Norma do M´aximo das Colunas na matriz em (3.2). Exemplo 3.2.1 Dado o sistema linear 8>: <> ¡7x1 +3x2 +2x3 = ¡2 x1 +3x2 x3 =3 ¡ x1 + x2 3x3 =¡¡1 vamos procurar uma aproxima¸c˜ao da solu¸c˜ao, com precis˜ao e =0:1 na Norma do M´aximo, usando o M´etodo de Gauss-Jacobi. Vamos tomar como condi¸c˜ao inicial o vetor x0 = (0:5, 0:5, 0:5)T . CAP´ITULO 3. SISTEMAS LINEARES 43 O primeiro passo ´e verificar se o crit´erio de convergˆencia ´e satisfeito. Calculando os ®k temos 15 ®1 =(a1;2+ a1;3)= < 1 ja1;1jjjjj7 12 ®2 =(a2;1+ a2;3)= < 1 ja 1 2;2jjjjj 23 ®3 =(a3;1+ a3;2)= < 1 ja3;3jjjjj3 Logo o crit´erio das linhas ´e satisfeito e com isto temos garantia que o M´etodo de Gauss-Jacobi converge. Montando o processo iterativo temos (k+1) 1 ¡2 - 3x2(k) - 2x(k) x = - 1 3 7 (k+1) 1 (k)(k) =+ x 2 ³3 - x1 3 (3.3) x 3 ³1¡- x (k+1) 1 (k)(k) x = - - x 3 1 2 3 1C Assimpara k =0 segueque Verificandoocrit´eriodeparadatemos 1C x(1) = ¡71 (¡2 - 3 * 0:5 - 2 * 0:5) = 0:642 1 x(1) = 1(3 - 0:5+0:5) = 1:000 (3.4) 2 3 x(1) = - 13 (¡1 - 0:5 - 0:5) = 0:666 3 . = . . jjx . B. . B. 0:642 - 0:500 1:000 - 0:500 0:142 10 10 0:500 =0:500 >e x - x = - x jj1 0:666 - 0:500 0:166 Para k = 1 temos x(2) 1 = - 1 7 (¡2 - 3 * 1:000 - 2 * 0:666) = 0:904 x(2) 2 = 1 3 (3 - 0:642 + 0:666) = 1:008 (3.5) x(2) = ¡31 (¡1 - 0:642 - 1) = 0:880 3 CAP´ITULO 3. SISTEMAS LINEARES 44 Verificando o crit´erio de parada temos . = . . jjx 1C 1C 0:262 = . B. . B. 0:904 - 0:642 1:008 - 1:000 21 21 0:008 =0:262 >e - x - x jj1 x 0:880 - 0:666 0:214 Devemos continuar as itera¸c˜oes at´e que o crit´erio de parada seja satisfeito (Exerc´icio). 3.2.4 M´etodo Iterativo de Gauss-Seidel A cada itera¸c˜ao xk se aproxima da solu¸c˜ao do sistema. Baseado nesta observa¸c˜ao vamos modificar o M´etodo de Gauss-Jacobi com o objetivo de acelerar a convergˆencia. Numa itera¸c˜ao k + 1, o M´etodo de Gauss-Jacobi calcula o elemento s pela equa¸c˜ao 0@ bs - s¡1 . j=1 as;jx (k) j - n j=s+1 as;jx (k) j 1A (3.6) 1 (k+1) = x s as;s k+1 k+1 k+1 Neste ponto os elementos x;x , ;x s¡1 , j´a foram calculados e espera-se que estes este 12 ¢¢· kk k jam mais pr´oximos da solu¸c˜ao que os elementos x1;x2, ;xAssim vamos substituir os ¢¢· s¡1. elementos da itera¸c˜ao k, que aparecem no primeiro somat´orio de (3.6), pelos correspondentes elementos da itera¸c˜ao k + 1, isto ´e 0@ bs - 1¡s. j=1 as;jx 1A n j=s+1 k+1 1 (k+1) (k) (k+1) = - x as;jx . j j s as;s Como estamos usando valores mais pr´oximos da solu¸c˜ao, o c´alculo de x ser´a mais preciso. s Enquanto >e fa¸jjjjx X 0@ Este procedimento ´e conhecido como M´etodo Iterativo de Gauss-Seidel, cujo o algoritmo ´e dado abaixo. Algoritmo: M´etodo Iterativo de Gauss-Seidel Input: Matriz A . IRn£n b, x0 . IRn e "> 0 k+1k - xca: Para s =1, 2;:::;n, fa¸ca: n s¡1 1. 1 (k+1) (k) (k+1) bs - - x as;jx as;jx j j s . as;s j=1 j=s+1 fim para fim enquanto Output: x . IRn : solu¸c˜ao do sistema CAP´ITULO 3. SISTEMAS LINEARES 45 3.2.5 Crit´erio de Sassenfeld O M´etodo de Gauss-Seidel tamb´em pode ser representado na forma matricial xk+1 = Cxk+g e crit´erios de convergˆencia podem ser obtidos impondo a condi¸c˜ao jjCj| < 1. Aplicando a Norma do M´aximo das Linhas obtemos o seguinte crit´erio de convergˆencia: Teorema 3.2.2 (Crit´erio de Sassenfeld) Dado o sistema linear Ax = b. Seja os ¯k de tal forma que: n 0. 1A k¡1 . j=1 1 ¯k = ¯j + < 1 para k =1, 2;:::;n jak;jj jak;jj jak;k| j=k+1 Ent˜ao o M´etodo de Gauss-Seidel gera uma seq¨uˆencia fxk} que converge para a solu¸c˜ao do sistema. Exemplo 3.2.2 Dado o sistema linear 8>: <> ¡7x1 +3x2 +2x3 = ¡2 x1 +2x2 x3 =2 ¡ x1 + x2 2x3 =0 - vamos procurar uma aproxima¸c˜ao da solu¸c˜ao, com precis˜ao e =0:1 na norma do m´aximo, usando o M´etodo de Gauss-Seidel. Vamos tomar como condi¸c˜ao inicial o vetor x0 = (0:5, 0:5, 0:5)T . O primeiro passo ´e verificar se o crit´erio de convergˆencia ´e satisfeito. Calculando os ¯k temos 15 ¯1 =(a1;2+ a1;3)= < 1 ja1;1jjjjj7 16 ¯2 =(a2;1¯1 + a2;3)= < 1 ja2;2jjjjj7 1 11 ¯3 =(a3;1¯1 + a3;2¯2)= < 1 ja3;3jjjjj14 Logo o crit´erio de Sassenfeld ´e satisfeito e com isto temos garantia que o M´etodo de Gauss- Seidel converge. Note que se aplicarmos o crit´erio das linhas obtemos que ®2 = ®3 =1, ou seja, o crit´erio das linhas n˜ao ´e satisfeito. Este ´e um exemplo em que n˜ao temos a garantia de convergˆencia do M´etodo de Gauss-Jacobi. Montando o processo iterativo temos (k+1) 1 ¡2 - 3x2(k) - 2x(k) x = - 1 3 7 (k+1) 1 (k+1) (k) =+ x 2 ³2 - x1 2 (3.7) x 3 x3(k+1) = ¡21 ³0 - x1(k+1) - x2(k+1) CAP´ITULO 3. SISTEMAS LINEARES 46 Assim para k =0 segue que x(1) = ¡71 (¡2 - 3 * 0:5 - 2 * 0:5) = 0:642 1 x(1) = 1(2 - 0:642 + 0:5) = 0:929 (3.8) 2 2 x(1) = - 1 (0 - 0:642 - 0:929) = 0:785 3 2 Verificando o crit´erio de parada temos = 1C )jj= xA. Para k =1 temos x(2) = 1(¡2 - 3 * 0:929 - 2 * 0:785) = 0:908 1 1C ¡7 x(2) = 1(2 - 0:908 + 0:785) = 0:938 (3.9) 2 2 x(2) = - 12 (0 - 0:908 - 0:938) = 0:923 3 Verificando o crit´erio de parada temos . B. 0:642 - 0:500 0:142 1 10 0:429 =0:429 >e x - x - x jj8 0 . B@ 0:929 - 0:500 0:785 - 0:500 0:285 1C 1C = . = . . jjx Devemos continuar as itera¸c˜oes at´e que o crit´erio de parada seja satisfeito (Exerc´icio). 3.3 Observa¸c˜oes Finais A escolha do m´etodo, que deve ser aplicado a um determinado problema, deve ser orientada nas caracter´isticas de cada m´etodo que apresentamos nesta se¸c˜ao. Os m´etodos diretos apresentam a solu¸c˜ao de qualquer sistema linear n˜ao singular, por´em n˜ao temos um controle sobre a precis˜ao da solu¸c˜ao. Aplicados em sistemas de grande porte e matriz cheia ( dimens˜ao acima 50 × 50 e poucos elementos ai;j = 0 ) apresentam grandes erros de arredondamentos. Os m´etodos iterativos permitem um controle sobre a precis˜ao da solu¸c˜ao, por´em estes n˜ao se aplicam a qualquer sistema. O sistema deve satisfazer certas condi¸c˜oes de convergˆencia para que determinado m´etodo seja aplicado. . B. 0:908 - 0:642 0:938 - 0:929 0:266 2 21 0:009 =0:266 >e x - x - x jj8 . B. 1 0:923 - 0:785 0:138 CAP´ITULO 3. SISTEMAS LINEARES 47 O M´etodo de Gauss-Jacobi ´e indicado para processamento paralelo ou vetorial, pois os elementos no passo k + 1 dependem somente dos elementos no passo k. Se T for o tempo que uma m´aquina seq¨uencial toma para executar uma itera¸c˜ao. Numa m´aquina paralela este tempo cai para dT=Npe, onde Np ´e o n´umero de processadores. O M´etodo de Gauss-Seidel n˜ao ´e indicado para processamento paralelo, pois o c´alculo de k+1 k+1 k+1 k+1 xs depende de x1 ;x2 , ;xs¡1 . Por´em este converge mais rapidamente que o M´etodo ¢¢· de Gauss-Jacobi, quando ambos s˜ao executado em processamento seq¨uencial. Al´em disso, todo sistema que satisfaz o Crit´erio das Linhas tamb´em satisfaz o Crit´erio de Sassenfeld. Ou seja, todo sistema em que podemos aplicar o M´etodo de Gauss-Jacobi, automaticamente podemos aplicar o M´etodo de Gauss-Seidel. 3.4 Exerc´icios Exerc´icio 3.1 Resolva o sistema linear abaixo, usando o M´etodo de Elimina¸c˜ao de Gauss. 8>: <> 1:7x1 +2:3x2 0:5x3 =4:55 1:1x1 +0:6x2 - 1:6x3 = ¡3:40 ¡ 2:7x1 0:8x2 +1:5x3 =5:50 - Exerc´icio 3.2 Ache a inversa da matriz abaixo. . B. 1 ¡22 2 ¡32 2 ¡21 . C. Exerc´icio 3.3 Um sistema ´e dito ser tridiagonal se este ´e formado pela diagonal principal, a primeira diagonal secund´aria inferior e a primeira diagonal secund´aria superior. Os outros elementos s˜ao nulos. Isto ´e, a matriz associada ´e da forma: 0. . a1;1 a1;2 0 0 0 · · · 0 0 0 a2;1 a2;2 a2;3 0 0 · · · 0 0 0 0 a3;2 a3;3 a3;4 0 · · · 0 0 0 0 0 a4;3 a4;4 a4;5 · · · 0 0 0 . . . . . . . . . . . . . . . . . . . . . . . . 0 0 0 0 0 0 0 0 0 0 · · · · · · an¡1;n¡2 0 an¡1;n¡1 an;n¡1 an¡1;n an;n 1. . Fa¸ca uma modifica¸c˜ao no M´etodo de Elimina¸c˜ao de Gauss explorando a forma do sistema. Implemente o algoritmo em MatLab. Exerc´icio 3.4 Considere o sistema linear . 0:0002x1 +2x2 =2 2x1 +2x2 =4 Calcule a solu¸c˜ao do sistema por Elimina¸c˜ao de Gauss e Pivotamento Parcial, usando 4 casas decimais, sem arredondamento. Calcule o res´iduo r = b - A¯x e comente seus resultados. CAP´ITULO 3. SISTEMAS LINEARES 48 Exerc´icio 3.5 Dado o sistema linear . 0:780x +0:563y =0:217 0:913x +0:659y =0:254 a) Calcule a solu¸c˜ao do sistema por (i)-Elimina¸c˜ao de Gauss e (ii)-Pivotamento Parcial, usando no m´inimo 7 casas decimais, sem arredondamento. b) Calcule o res´iduo r = b - A¯x para os casos (i) e (ii). c) Se no ´item a) tiv´essemos usado 3 casas decimais, o que ocorreria com a solu¸c˜ao do sistema? Comente seus resultados. Exerc´icio 3.6 Mostre que, se um sistema linear satisfaz o Crit´erio das Linhas, ent˜ao este tamb´em satisfaz o Crit´erio de Sassenfeld. Exerc´icio 3.7 Seja k um n´umero inteiro, positivo, considere: 8.. .. kx1 + x2 =2 k kx1 +2x2 + x3 =3 5 kx1 + x2 +2x3 =2 a) Verifique para que valores de k , a convergˆencia do M´etodo de Gauss-Jacobi pode ser garantida. b) Verifique para que valores de k, a convergˆencia do M´etodo de Gauss-Seidel pode ser garantida. c) Utilize um m´etodo iterativo adequado para calcular a aproxima¸c˜ao da solu¸c˜ao deste sistema de equa¸c˜oes considerando: (i) x(0) = (1:0, 1:0, 1:0)T (ii) Escolha k como o menor inteiro que satisfa¸ca as condi¸c˜oes de convergˆencia. (iii) Fa¸ca duas itera¸c˜oes e calcule o erro absoluto cometido, usando a norma do m´aximo. Exerc´icio 3.8 Dado o sistema Ax = b podemos montar um processo iterativo da seguinte forma x k+1 =(I + A)x k - b a) Enuncie uma condi¸c˜ao suficiente de convergˆencia baseada na Norma do M´aximo das Linhas. b) Fa¸ca trˆes itera¸c˜oes do processo acima para o sistema linear ¡1:3x1 +0:3x2 =1 (0) . 0:8 tomando x= 0:5x1 0:5x2 =0 0:8 - Cap´itulo 4 Ajuste de Curvas: M´etodo dos M´inimos Quadrados Em geral, experimentos em laborat´orio gera uma gama de dados que devem ser analisados para a cria¸c˜ao de um modelo. Obter uma fun¸c˜ao matem´atica que represente (ou que ajuste) os dados permite fazer simula¸c˜oes do processo de forma confi´avel, reduzindo assim repeti¸c˜oes de experimentos que podem ter um custo alto. Neste cap´itulo vamos analisar o esquema dos M´inimos Quadrados, que fornece uma fun¸c˜ao que melhor represente os dados. 4.1 M´etodo dos M´inimos Quadrados -Caso Discreto Dado um conjunto de pontos (xk;f(xk)), k =0, 1, 2, :::, m. O problema de ajuste de curvas consiste em encontrar uma fun¸c˜ao '(x) tal que o desvio em cada ponto k, definido por dk = f(xk) - '(xk), seja m´inimo, onde '(x) ´e uma combina¸c˜ao linear de fun¸c˜oes cont´inuas gi(x);i =1, 2, :::, n, escolhidas de acordo com os dados do problema. Isto ´e '(x)= ®1g1(x)+ ®2g2(x)+ + ®ngn(x) ¢¢· Neste caso dizemos que o ajuste ´e linear. A escolha das fun¸c˜oes gi depende do gr´afico dos pontos, chamado de diagrama de dispers˜ao, atrav´es do qual podemos visualizar que tipo de curva que melhor se ajusta aos dados. Exemplo 4.1.1 Vamos considerar a tabela de pontos dada abaixo. x 0:10 0:20 0:50 0:65 0:70 0:80 0:90 1:10 1:23 1:35 1:57 1:70 1:75 1:80 1:94 f(x) 0:19 0:36 0:75 0:87 0:91 0:96 0:99 0:99 0:94 0:87 0:67 0:51 0:43 0:36 A an´alise do gr´afico de dispers˜ao (Fig. 4.1) mostra que a fun¸c˜ao que procuramos se comporta como uma par´abola. Logo poder´iamos escolher as fun¸c˜oes g1(x)=1;g2(x)= x e 49 ´ CAP´ITULO 4. AJUSTE DE CURVAS: M ETODO DOS M ´INIMOS QUADRADOS 50 Figura 4.1: Diagrama de Dispers˜ao g3(x)= x2, pois '(x)= ®1g1(x)+ ®2g2(x)+ ®3g3(x) representa “todas” as par´abolas e com a escolha adequada dos ®i teremos aquela que melhor se ajusta aos pontos. O M´etodo dos M´inimos Quadrados consiste em determinar os ®i de tal forma que a soma dos quadrados dos desvios em seja m´inimo, Isto ´e: Achar os ®i que minimizam a fun¸c˜ao '(xk) m F (®1;®2;:::;®n)= . [f(xk)- z(®1g1(xk)+}. ®ngn(xk))]. 2 . k=1 ¢¢· A fun¸c˜ao F ´e uma fun¸c˜ao que satisfaz F (®) = 0 8a . IRm . Isto ´e, uma fun¸c˜ao limitada inferiormente e portanto esta tem um ponto de m´inimo. Este ponto pode ser determinado pelo teste da primeira derivada, sendo @F . =0 i =1, . . . , n. @®i (®1;:::;®n) Desta forma temos m . (xk)] gj(xk)=0 ¡2 k=1 [f(xk) - ®1g1(xk) - ®2g2(xk) ¡¢¢· ®ngn. ´ CAP´ITULO 4. AJUSTE DE CURVAS: M ETODO DOS M ´INIMOS QUADRADOS 51 8. ®1 mmm . k=1 k=1 k=1 k=1 mmm. m g1(xk)g1(xk)+ ®2 g1(xk)g2(xk)+ g1(xk)gn(xk) f(xk)g1(xk) + ®n = ¢¢· m g2(xk)g1(xk)+ ®2 <. g2(xk)g2(xk)+ ¢¢· + ®n g2(xk)gn(xk) = f(xk)g2(xk) ®1 k=1 k=1 k=1 k=1 ..... . ..... . ..... . mmm. k=1 k=1 k=1 k=1 (xk)g1(xk)+ ®2 m (xk)g2(xk)+ (xk)gn(xk) f(xk)gn(xk) + ®n ®1 = gn gn gn : ¢¢· Que representa um sistema linear n × n da forma 8. <. . a1;1®1 + a1;2®2 + a1;3®3 a1;n®n = b1 ¢¢· a2;1®1 + a2;2®2 + a2;3®3 a2;n®n = b2 ¢¢· a3;1®1 + a3;2®2 + a3;3®3 a3;n®n = b3 ¢¢· ... .. ... .. ... .. an;1®1 + an;2®2 + an;3®3 an;n®n = bn ¢¢· onde mm. k=1 k=1 Este sistema tem uma ´unica solu¸c˜ao se os vetores formados por gk =(gk(x1);gk(xn)) ¢¢· s˜ao linearmente independentes. Isto ´e equivalente a ter as fun¸c˜oes gi(x) linearmente inde pendentes. A matriz A associada ao sistema ´e uma matriz sim´etrica, ou seja ai;j = aj;i. Logo, para um sistema n × n, ser´a necess´ario calcular (n2 - n)=2 elementos. Exemplo 4.1.2 Usando a tabela do exemplo (4.1.1), vamos ajustar os dados por uma par´abola. Para isto vamos tomar g1(x)=1;g2(x)= x e g3(x)= x2 . Calculando cada uma das fun¸c˜oes nos pontos xk temos. gi(xk)gj(xk)e bi f(xk)gi(xk) ai;j = = x 0:10 0:20 0:50 0:65 0:70 0:80 0:90 1:10 1:23 1:35 1:57 1:70 1:75 1:80 1:94 f(x) 0:19 0:36 0:75 0:87 0:91 0:96 0:99 0:99 0:94 0:87 0:67 0:51 0:43 0:36 0:11 g1(x) 1:0 1:0 1:0 1:0 1:0 1:0 1:0 1:0 1:0 1:0 1:0 1:0 1:0 1:0 1:0 g2(x) 0:10 0:20 0:50 0:65 0:70 0:80 0:90 1:10 1:23 1:35 1:57 1:70 1:75 1:80 1:94 g3(x) 0:01 0:04 0:25 0:42 0:49 0:64 0:81 1:21 1:51 1:82 2:46 2:89 3:06 3:24 3:76 15. Calculando os elementos da matriz e o vetor dos termos independentes temos g1(xk) * g1(xk) = 15 a1;1 = k=1 15X g1(xk) * g2(xk) = 16:29 = a2;1 a1;2 = k=1 15X g1(xk) * g3(xk) = 22:62 = a3;1 a1;3 = k=1 ´ CAP´ITULO 4. AJUSTE DE CURVAS: M ETODO DOS M ´INIMOS QUADRADOS 52 15 a2;2 = . g2(xk) * g2(xk) = 22:62 k=1 15 a2;3 = . g2(xk) * g3(xk) = 34:92 = a3;2 k=1 15 a3;3 = . g3(xk) * g3(xk) = 57:09 k=1 15 b1 = . f(xk) * g1(xk)=9:91 k=1 15 b2 = . f(xk) * g2(xk) = 10:28 k=1 15 b3 = . f(xk) * g3(xk) = 12:66 k=1 Obtendo assim um sistema linear que pode ser resolvido por um esquema num´erico estudado no cap´itulo 3. A solu¸c˜ao do sistema ´e dado por ®1 =0:00;®2 =1:99;®3 = ¡0:99 Portanto a fun¸c˜ao . ´e dada por '(x)=1:99x - 0:99x 2 A figura 4.2 compara a fun¸c˜ao '(x) com o gr´aficos dos pontos. Atrav´es da fun¸c˜ao . podemos aproximar valores de f em pontos que n˜ao pertencem a tabela. Por exemplo: f(0) '(0) = 0 ˜ f(1) '(1) = 1 ˜ f(2) '(2) = 0:02 ˜ No exemplo ajustamos os dados a uma par´abola, mas outras fun¸c˜oes bases poderiam ser usadas. Como exemplo, poder´iamos pensar que os dados representam o primeiro meio ciclo de uma fun¸c˜ao senoidal. E neste caso poder´iamos tomar g1(x)=1e g2(x) = sen(2¼x) Mas afinal qual seria a melhor escolha? (Veja exerc´icio 4.1) O desvio fornece uma medida que pode ser usada como parˆametro de compara¸c˜ao entre ajustes diferentes. No caso do ajuste pela par´abola temos que o desvio ´e dado por 15 D = X(f(xk) - '(xk))2 =0:0019 k=1 Se o ajuste feito por uma fun¸c˜ao senoidal tiver um desvio menor, ent˜ao este ajuste representaria melhor os dados. Outro ponto a ser observado ´e que a dimens˜ao do sistema linear depende do n´umero de fun¸c˜oes bases que estamos usando. No caso da par´abola usamos trˆes fun¸c˜oes bases e temos um sistema 3£3. No caso de uma fun¸c˜ao senoidal teremos um sistema 2 × 2. ´ CAP´ITULO 4. AJUSTE DE CURVAS: M ETODO DOS M ´INIMOS QUADRADOS 53 Figura 4.2: Diagrama de Dispers˜ao com o gr´afico da '(x). 4.2 M´etodo dos M´inimos Quadrados -Caso Cont´inuo No caso cont´inuo temos uma fun¸c˜ao f(x) dada num intervalo [a, b] e n˜ao mais uma tabela de pontos. O procedimento ´e an´alogo ao caso discreto. Escolhidas as fun¸c˜oes bases gi devemos determinar a fun¸c˜ao '(x)= ®1g1(x)+ ®2g2(x)+ + ®ngn(x) de modo que o desvio seja ¢¢· m´inimo, onde D = . b (f(x) - '(x))2dx a Neste caso os ®i tamb´em s˜ao determinados pela resolu¸c˜ao de um sistema, onde o elemento ai;j e ´e obtido atrav´es do produto interno entre as fun¸c˜oes gi(x)e gj(x) e o elementos bi pelo produto interno entre f(x)e gi(x), sendo . b . b ai;j = gi(x)gj(x)dx e bi = f(x)gi(x)dx aa Exemplo 4.2.1 Vamos determinar a melhor par´abola que se ajuste a fun¸c˜ao f(x)= sen(¼x) no intervalo [0, 1]. Para isto devemos tomar, como fun¸c˜oes bases, as fun¸c˜oes g1(x)=1, g2(x)= x e g3(x)= x2 . Calculando os termos do sistema linear temos . 1 a1;1 = g1(x)g1(x)dx =1 0 . 1 1 a1;2 = g1(x)g2(x)dx == a2;1 0 2 ´ CAP´ITULO 4. AJUSTE DE CURVAS: M ETODO DOS M ´INIMOS QUADRADOS 54 . 1 1 a1;3 = g1(x)g3(x)dx == a3;1 0 3 . 1 1 a2;2 = g2(x)g2(x)dx = 0 3 . 1 1 a2;3 = g2(x)g3(x)dx == a3;2 0 4 . 1 1 a3;3 = g3(x)g3(x)dx = 0 25 . 1 b1 = f(x)(x)g1(x)dx =0:636 0 . 1 b2 = f(x)g2(x)dx =0:318 0 . 1 b3 = f(x)g3(x)dx =0:189 0 cuja a solu¸c˜ao ´e dada por ®1 = ¡0:027, ®2 =4:032 e ®3 = ¡4:050. Assim temos que '(x)= ¡0:027 + 4:032x - 4:050x2 . A figura (4.3) mostra o gr´afico comparativo entre a fun¸c˜ao f(x) ( linha: —) e o ajuste '(x) (linha: ). ¢¢¢ Figura 4.3: — : f(x) = sen(¼x); : '(x) ¢¢· ´ CAP´ITULO 4. AJUSTE DE CURVAS: M ETODO DOS M ´INIMOS QUADRADOS 55 4.3 Ajuste N˜ao Linear Existem casos, onde o diagrama de dispers˜ao de uma fun¸c˜ao indica que os dados devem ser ajustado por uma fun¸c˜ao que n˜ao ´e linear com rela¸c˜ao aos ®i. Como exemplo, vamos considerar os dados 0 0:5 1:0 1:5 2:0 2:5 3:0 ¡1:0 ¡0:5 f(x) 0:157 0:234 0:350 0:522 0:778 1:162 1:733 2:586 3:858 Montando o diagrama de dispers˜ao (Veja figura 4.4) podemos considerar que f(x) tem um comportamento exponencial. Isto ´e, f(x) ˜ '(x)= ®1e®2x . Desta forma, um processo Figura 4.4: Diagrama de dispers˜ao de lineariza¸c˜ao deve ser empregado, para que seja poss´ivel aplicar o M´etodo dos M´inimos Quadrados. Neste caso podemos proceder da seguinte forma. f(x)= ®1e ®2x z = ln(f(x)) = ln(®1)+ ®2x. . Fazendo ¯1 = ln(®1)e ¯2 = ®2 o problema consiste em ajustar os dados de z por uma reta. Para isto tomamos g1(x)=1e g2(x)= x. Calculando as fun¸c˜oes em cada um dos pontos ´ CAP´ITULO 4. AJUSTE DE CURVAS: M ETODO DOS M ´INIMOS QUADRADOS 56 temos x ¡1:0 ¡0:5 0 0:5 1:0 1:5 2:0 2:5 3:0 f(x) 0:157 0:234 0:350 0:522 0:778 1:162 1:733 2:586 3:858 z = ln(f(x)) ¡1:851 ¡1:452 ¡1:049 ¡0:650 ¡0:251 0:150 0:549 0:950 1:350 g1(x) 1:000 1:000 1:000 1:000 1:000 1:000 1:000 1:000 1:000 g2(x) 0 0:5 1:0 1:5 2:0 2:5 3:0 ¡1:0 ¡0:5 Calculando os elementos da matriz e vetor dos termos independente temos que 9 a1;1 = . g1(xk) * g1(xk)=9 k=1 9 a1;2 = . g1(xk) * g2(xk)=9= a2;1 k=1 9 a2;2 = . g2(xk) * g2(xk) = 24 k=1 15 b1 = . z(xk) * g1(xk)= ¡2:254 k=1 15 b2 = . z(xk) * g2(xk)=9:749 k=1 Cuja a solu¸c˜ao ´e dada por ¯1 = ¡1:050 e ¯2 =0:800 Desta forma os valores de ®i s˜ao dados por: ®1 = e ¯1 =0:349 e ®2 = ¯2 =0:800 Portanto temos '(x)= ®1e ®2 =0:349e 0:800x A figura 4.5 mostra a compara¸c˜ao dos dados com a fun¸c˜ao obtida. Para verificar se fun¸c˜ao escolhida para a aproxima¸c˜ao foi bem feita, usamos o teste de alinhamento. Este consiste em tomarmos os dados “linearizados”, isto ´e, os pontos z da tabela, e fazer o diagrama de dispers˜ao. Se os pontos estiverem alinhados, ent˜ao a escolha da fun¸c˜ao foi boa. A figura 4.6 mostra o diagrama de dispers˜ao dos dados em z, obtidos no nosso exemplo. Podemos concluir que a nossa escolha pela exponencial foi uma escolha acertada. 4.4 Observa¸c˜oes Finais O ajuste de curvas permite extrapolar os valores tabelados. Isto ´e, se os dados est˜ao tabelados num intervalo [x0;xm] podemos aproximar um ¯x 6 . [x0;xm] com uma certa seguran¸ca. Como ´ CAP´ITULO 4. AJUSTE DE CURVAS: M ETODO DOS M ´INIMOS QUADRADOS 57 Figura 4.5: Diagrama de Dispers˜ao e o Gr´afico da '(x) os dados prov´em de experimentos que est˜ao sujeitos a erros de medi¸c˜oes, podemos ter mais de um valor para um determinado ponto. (Veja exerc´icio 4.4) A fun¸c˜ao obtida considera os dois valores faz “ uma m´edia ” entre estes valores. Os elementos ai;j s˜ao obtidos pelo produto interno entre as fun¸c˜oes gi e gj definidos por m Caso Discreto: hgi;gj. = . gi(xk)gj(xk) k=1 . b Caso Cont´inuo: hgi;gj. = gi(x)gj(x)dx a Se as fun¸c˜oes gi forem ortogonais, isto ´e . 0, para i = j hgi;gj. = 6 ki, para i = j a matriz obtida ser´a diagonal e conseq¨uentemente a solu¸c˜ao do sistema ´e imediata. Como exemplo, veja exerc´icio 4.5. ´ CAP´ITULO 4. AJUSTE DE CURVAS: M ETODO DOS M ´INIMOS QUADRADOS 58 Figura 4.6: Diagrama dos dados linearizados 4.5 Exerc´icios Exerc´icio 4.1 Usando os dados abaixo, fa¸ca um ajuste de curva com g1(x)=1 e g2(x)= sen(2¼x). Calcule o desvio e compare com os resultados obtidos no Exemplo (4.1.1). x 0:10 0:20 0:50 0:65 0:70 0:80 0:90 1:10 1:23 1:35 1:57 1:70 1:75 1:80 1:94 f(x) 0:19 0:36 0:75 0:87 0:91 0:96 0:99 0:99 0:94 0:87 0:67 0:51 0:43 0:36 0:11 Exerc´icio 4.2 Dada a tabela abaixo x 0.50 0.75 1.00 1.50 2.00 2.50 3.00 f(x) 0.479 0.681 0.841 0.997 0.909 0.598 0.141 Entre os grupos de fun¸c˜oes bases abaixo, escolha aquele que representar´a melhor re sultado num ajuste de curvas. Justifique sua escolha. Fa¸ca um ajuste considerando sua escolha. Grupo I: g1(x)=1 e g2(x)= x. Grupo II: g1(x)=1 e g2(x)= ex . Grupo III: g1(x)=1 e g2(x)= x2 . ´ CAP´ITULO 4. AJUSTE DE CURVAS: M ETODO DOS M ´INIMOS QUADRADOS 59 Exerc´icio 4.3 A tabela abaixo representa o calor espec´ifico da ´agua em fun¸c˜ao da temperatura. t(oC) 0 5 10253035 C(t) 1.00762 1.00392 1.00153 0.99852 0.99826 0.99818 Fac ca um ajuste linear, um quadr´atico e um c´ubico. Fa¸ca um ajuste n˜ao linear da forma '(x)= ®1e¡®2 , com ®1;®2 > 0. Calcule o desvio e ache uma aproxima¸c˜ao para t = 15 em cada um dos casos. Sabendo que o valor exato da fun¸c˜ao C(15)= 1:00000, qual dos casos acima apresentou melhor aproxima¸c˜ao? 1 Exerc´icio 4.4 Ajustar os dados abaixo `a fun¸c˜ao z = 1+ e(®1x+®2) x 0.1 0.3 0.5 0.5 0.7 0.8 0.8 1.10 1.30 1.80 f(x) 0.833 0.625 0.500 0.510 0.416 0.384 0.395 0.312 0.277 0.217 Verifique, pelo teste do alinhamento, qual ´e a melhor escolha para ajustar os dados entre as fun¸c˜oes z = ®1®2 x e z = ®1x®2 .(Obs: Note que neste caso a tabela apresenta dois valores diferentes para os pontos x =0:5 e x =0:8. Isto pode ocorrer se estamos considerando que os dados prov´em de experimentos que s˜ao realizados mais de uma vez. ) 1 Exerc´icio 4.5 Usando os polinˆomios de Legendre g1(x)=1, g2(x)= x e g2(x)= (3x 2 ¡1), 2 que s˜ao ortogonais em [¡1, 1], ache a melhor par´abola que aproxima a fun¸c˜ao f(x) = cos(x) no intervalo [¡3, 3].(Obs: A ortogonalidade dos polinˆomios nos forneceria uma matriz diagonal, se o ajuste fosse feito no intervalo [¡1, 1]. Logo devemos fazer uma mudan¸ca de vari´avel de tal forma que obteremos novas gi que ser˜ao ortogonais em [¡3, 3] ) Cap´itulo 5 Interpola¸c˜ao Polinomial A interpola¸c˜ao ´e outra forma de encontrar uma fun¸c˜ao que represente um conjunto de dados tabelados. Interpolar um conjunto de dados (xk;fk), k =0, 1, ;n, consiste em encontrar ¢¢· uma fun¸c˜ao pn(x), escolhida numa classe de fun¸c˜oes, tal que esta satisfa¸ca certas propriedades. Neste cap´itulo vamos considerar o caso onde pn(x) ´e um polinˆomio de tal forma que fk = p(xk);k =0, 1, 2, ;n. ¢¢· Esta condi¸c˜ao ´e chamada de condi¸c˜ao de interpola¸c˜ao e o polinˆomio que satisfaz esta condi¸c˜ao ´e chamado de polinˆomio interpolador. Teorema 5.0.1 (Existˆencia e Unicidade) Dado o conjunto de n +1 pontos distintos ¢¢· 66unico polinˆ (xk;fk);k =0, 1;n, Isto ´e, xk = xj para k = j. Existe um ´omio p(x) de grau menor ou igual a n, tal que p(xk)= fk para k =0, 1, 2, ;n. ¢¢· n Prova: Seja p(x)= a0 + a1x + a2x2 ++ anx. Para obter os ai usamos a condi¸c˜ao de ¢¢· interpola¸c˜ao fk = p(xk) para k =0, 1, 2, ;n. Logo, segue que: ¢¢· f0 = p(x0)= a0 + a1x0 + a2x02 ++ anx0 n ¢¢· f1 = p(x1)= a0 + a1x1 + a2x12 ++ anx1 n ¢¢· . ... . ... .=.. . fn = p(xn)= a0 + a1xn + a2xn 2 ++ anxn n ¢¢· Que corresponde ao sistema linear da forma 0. 1 x0 x02 x0 n ¢¢· 1 x1 x12 x1 n ¢¢· 0. 1. a0 a1 0. 1. f0 f1 1. 1 x2 x22 x2 n a2 = f2 ¢¢· ..... . ..... . @ ... . @A . @A . A 1 xn xn 2 xn n ¢¢· an fn 60 ˜ CAP´ITULO 5. INTERPOLAC¸ AO POLINOMIAL 61 A matriz A, associada ao sistema, ´e uma matriz de Vandermonde, cujo o determinante ´e dado por nl¡1 Det(A)= (xl - xj) l=1 j=0 Como xl = xj para l = j, segue que o determinante da matriz A ´e diferente de zero e 66unica solu¸aoportanto o sistema admite uma ´c˜ Exemplo 5.0.1 Vamos achar uma aproxima¸c˜ao para f(0:3) usando o polinˆomio interpolador dos dados abaixo. 0.0 0.2 0.4 xk 4.00 3.84 3.76 fk Como temos, trˆes pontos (n +1 = 3), o grau do polinˆomio ser´a menor ou igual a dois. Logo p(x)= a0 + a1x + a2x 2 Impondo a condi¸c˜ao fk = p(xk) obtemos: f0 =4:00 = p(0) = a0 + a10+ a202 f1 =3:84 = p(0:2) = a0 + a10:2+ a20:22 f2 =3:76 = p(0:4) = a0 + a10:4+ a20:42 Que equivale ao sistema linear na forma matricial . B. 10 0 4:00 10:20:04 3:84 10:40:16 3:76 . C. A solu¸c˜ao deste sistema ´e a0 =4, a1 = ¡1 e a2 =1, obtendo assim p(x)= x 2 - x +4 Desta forma f(0:3) ˜ p(0:3) = 3:79 Existem outras formas de encontrar o polinˆomio interpolador que a resolu¸c˜ao de sistemas. Teoricamente estes procedimentos resultam no mesmo polinˆomio pn(x). A escolha de uma ou outra forma depende dos dados que devemos interpolar. ˜ CAP´ITULO 5. INTERPOLAC¸ AO POLINOMIAL 62 5.1 Forma de Lagrange Vamos considerar o conjunto de n+1 pontos (xk;fk), k =0, 1;:::n distintos e vamos considerar o polinˆomio representado por n pn(x)= f0L0(x)+ f1L1(x)+ + fnLn(x)= . fkLk(x) ¢¢· k=0 onde Lk(x) ´e um polinˆomio de grau = n que satisfaz a rela¸c˜ao . 0 se k = j Lk(xj)= 6 1 se k = j Com isto temos que 0 0 + ¢¢· 1 %fL ()+ xjjj + ¢¢· 0 + fnLn%(xj)= fj pn (xj)= f0L0%(xj)+ f1L1%(xj) Logo pn(x) ´e o polinˆomio interpolador de f(x) nos pontos x0;x1;:::;xn. Os polinˆomios Lk(x) s˜ao chamados de polinˆomios de Lagrange e estes s˜ao obtidos da forma: Lk(x)= (x - x0)(x - x1)(x - xk¡1)(x - xk+1)(x - xn) ¢¢· ¢¢· (xk - x0)(xk - x1)(xk - xk¡1)(xk - xk+1)(xk - xn) ¢¢· ¢¢· Exemplo 5.1.1 Vamos considerar a tabela de pontos do exemplo anterior x 0.0 0.2 0.4 f(x) 4.00 3.84 3.76 Calculando os Lk(x) temos (x - x1)(x - x2)(x - 0:2)(x - 0:4) 1 L0(x)= (x0 - x1)(x0 - x2) = (0 - 0:2)(0 - 0:4) = 0:08(x 2 - 0:6x +0:08) = (x - x0)(x - x2) = (x - 0)(x - 0:4) = ¡1(x 2 L1(x) (x1 - x0)(x1 - x2) (0:2 - 0)(0:2 - 0:4) 0:04- 0:4x) (x - x0)(x - x1)(x - 0)(x - 0:2) 1 L2(x)= (x2 - x0)(x2 - x1) = (0:4 - 0)(0:4 - 0:2) = 0:08(x 2 - 2:6x) Assim temos que p(x)= x 2 - x +4 Observe que o polinˆomio ´e o mesmo obtido pela resolu¸c˜ao de sistema. Isto j´a era esperado, pois o polinˆomio interpolador ´´ e unico. ˜ CAP´ITULO 5. INTERPOLAC¸ AO POLINOMIAL 63 5.2 Forma de Newton A forma de Newton do polinˆomio interpolador ´e baseada nos operadores de diferen¸cas divididas. Seja f(x) uma fun¸c˜ao tabelada em n + 1 pontos distintos x0;x1;:::;xn. Definimos o operador de diferen¸ca dividida de ordem zero em xk por: f[xk]= f(xk). O operador de diferen¸ca dividida de ordem um, nos pontos xk;xk+1, ´e definido da seguinte forma: f[xk;xk+1]= f[xk] - f[xk+1] xk - xk+1 Este valor pode ser interpretado como uma aproxima¸c˜ao para a primeira derivada de f(x), em xk. O operador de diferen¸ca dividida de ordem dois, nos pontos xk;xk+1;xk+2, ´e definido da seguinte forma: f[xk;xk+1;xk+2]= f[xk;xk+1] - f[xk+1;xk+2] . xk - xk+2 De forma an´aloga, definimos o operador diferen¸ca dividida de ordem n, nos pontos xk;xk+1;:::;xk+n, da seguinte forma: f[xk;xk+1;:::;xk+n]= f[xk;xk+n¡1] - f[xk+1;xk+n] . xk - xk+n Note que a forma de c´alculo desses operadores ´e construtiva, no sentido de que para obter a diferen¸ca dividida de ordem n necessitamos das diferen¸cas divididas de ordem n - 1;n ¡2;:::, 1, 0. Um esquema pr´atico para o c´alculo desses operadores ´e dado pela tabela abaixo x x0 x1 x2 f[xk] f[xk;xk+1] f[xk;xk+1;xk+2] f[xk;xk+1;:::;xk+n ¢¢· ] f0 > f0¡f1 f1 x0¡x1 > f[x0;x1]¡f[x1;x2] f1¡f2 x0¡x2 > f2 x1¡x2 > f[x1;x2]¡f[x2;x3] x1¡x3 > f2¡f3 x2¡x3 > f[x0;:::;xn¡1]¡f[x1;:::;xn] x0¡xn x3 f3 .. . .. . .. . xn¡1 fn¡1 > fn¡1¡fn xn¡1¡xn xn fn ˜ CAP´ITULO 5. INTERPOLAC¸ AO POLINOMIAL 64 5.2.1 Constru¸c˜ao do Polinˆomio Vamos considerar o conjunto de pontos x0;x1;:::;xn, onde conhecemos os valores da fun¸c˜ao f(x), dados por f0;f1;:::;fn. Calculando a diferen¸ca dividida de ordem dois entre os pontos x, x0;x1 temos f[x, x0;x1]= f[x, x0] - f[x0 - x1] x - x1 Isolando a diferen¸ca de ordem um que depende de x segue que f[x, x0]= f[x0;x1]+(x - x1)f[x, x0;x1] Aplicamos a defini¸c˜ao de diferen¸ca de ordem um no primeiro termo, assim segue que f(x) - f(x0) = f[x0;x1]+(x - x1)f[x, x0;x1], x - x0 e isto implica que f(x)= f(x0)+ x - x0f[x0;x1]+(x - x0)(x - x1)f[x, x0;x1]= p1(x)+ E1(x). Ou seja a fun¸c˜ao f(x) ´e igual a um polinˆomio de grau um ,p1(x), mais uma fun¸c˜ao E1(x) que depende da diferen¸ca dividida de ordem dois. Desta forma podemos dizer que a fun¸c˜ao f(x) ´e aproximada por p1(x) com erro de E1(x). O polinˆomio p1(x) ´e o polinˆomio interpolador de f(x), nos pontos x0;x1, pois, 0 0 p1(x0)= f(x0)+(x0 ¡%x0)f[x0;x1]+(x0 ¡%x0)(x - x1)f[x, x0;x1] = f(x0) 0 p1(x1)= f(x0)+(x1 - x0)f[x0;x1]+(x1 - x0)(x ¡%x1)f[x, x0;x1] = f(x0)+(x1 - x0) f(x0) - f(x1) x0 - x1 = f(x1) De forma an´aloga, podemos calcular a diferen¸ca dividida de ordem n, sobre os pontos x, x0;x1;:::;xn, obtendo f(x)= pn(x)+ En(x), onde pn(x)= f(x0)+(x - x0)f[x0;x1]+(x - x0)(x - x1)f[x0;x1;x2]+ + ¢¢· (x - x0)(x - x1)(x - xn¡1)f[x0;x1 :::;xn] (5.1) ¢¢· En(x)=(x - x0)(x - x1)(x - xn)f[x, x0;x1;:::;xn] (5.2) ¢¢· ˜ CAP´ITULO 5. INTERPOLAC¸ AO POLINOMIAL 65 Assim podemos aproximar f(x) por pn(x), sendo que o erro ´e dado por En(x). O polinˆomio pn(x) ´e o polinˆomio interpolador de f(x) sobre os pontos x0;x1;:::xn, pois p(xj)= f(xj), para j =0, 1;:::;n. Exemplo 5.2.1 Vamos considerar a fun¸c˜ao f(x) tabelada abaixo. x 00:511:5 f(x) 0:0000 1:1487 2:7183 4:9811 Montando a tabela das diferen¸cas divididas temos xkf[xk]f[xk;xk+1]f[xk;::;xk+2]f[xk;::;xk+3] 0:00:00002:29740:51:14870:84183:13920:363061:02:71831:38644:52561:54:9811 Atrav´es da equa¸c˜ao (5.1) podemos notar que as diferen¸cas divididas que necessitamos s˜ao as primeiras de cada coluna. Logo polinˆomio interpolador ´e dado por p(x)= f(x0)+ x - x0f[x0;x1]+(x - x0)(x - x1)f[x0;x1;x2]+(x - x0)(x - x1)(x - x2)f[x0;x1;x2;x3] =0:00 + (x - 0:0)2:2974 + (x - 0:0)(x - 0:5)0:8418 + (x - 0:0)(x - 0:5)(x - 1:0)0:36306 =2:05803x +0:29721x 2 +0:36306x 3 5.3 Estudo do Erro A equa¸c˜ao (5.2) representa o erro cometido na interpola¸c˜ao sobre os pontos x0;:::;xn. Se aproximamos f(¯x) ˜ pn(¯x) o erro ser´a dado por En(¯x). Por´em este depende da diferen¸ca dividida f[¯x, x0;x1;:::;xn], que por sua vez, depende do valor de f(¯x). Como a fun¸c˜ao f(x) ´e tabelada, n˜ao temos como calcular este valor. Estimativas para o erro podem ser obtidas se conhecemos algumas propriedades da fun¸c˜ao. Teorema 5.3.1 Sejam x0;x1;:::;xn, n +1 pontos distintos. Seja pn(x) o polinˆomio interpolador sobre x0;x1;:::;xn. Se f(x) ´e n +1 vezes diferenci´avel em [x0;xn] ent˜ao para x¯. [x0;xn] o erro ´e dado por: f(n+1)(») En(¯x)= f(¯x) - pn(¯x) = (¯x - x0)(¯x - x1) ¢¢· (¯x - xn) (n + 1)! com . . [x0;xn] ˜ CAP´ITULO 5. INTERPOLAC¸ AO POLINOMIAL 66 Prova: Seja G(x)=(x - x0)(x - x1)(x - xn), logo G(xi) = 0 para i =0, 1;:::n. ¢¢· Seja H(t)= En(x)G(t) - En(t)G(x), logo H satisfaz 1: H(t) possui derivadas at´e ordem n + 1, pois G e En possuem derivadas at´e esta ordem. 2: H(t) possui pelo menos (n + 2) zeros em [x0;xn], pois para t = xi temos 00 H(xi)= En(x)G(xi)- En(xi)G(x)=0 e para t = x temos H(x)= En(x)G(x) - En(x)G(x) = 0. 3: Aplicando o Teorema de Rolle a H(t) e suas derivadas at´e ordem n + 1, temos H(t) tem n + 2 zeros em [x0, xn] H0(t) tem n + 1 zeros em [x0, xn] H00(t) tem n zeros em [x0, xn] . . . . . . H(n+1) tem 1 zeros em [x0, xn] Por ouro lado temos que H(n+1)(t)= En(x)G(n+1)(t) - En (n+1)(t)G(x) onde, En (n+1)(t)= f(n+1)(t) - pn (n+1)(t) G(n+1)(t)=(n + 1)! Como o polinˆomio pn ´e de grau n temos que pn (n+1)(t) = 0 e segue que H(n+1)(t)= En(x)(n + 1)! - f(n+1)(t)G(x) A fun¸c˜ao H(n+1)(t) possui um zero em [x0;xn] que vamos chamar de ». Substituindo na equa¸c˜ao acima temos que f(n+1)(») En(x)=(x - x0)(x - x1) ¢¢· (x - xn) (n + 1)! Na pr´atica usamos um limitante para o erro, sendo En(x)(x - x0)(x - x1)(x - xn)max jf(n+1)(x)j;jj·j¢¢· | x2[x0;xn] (n + 1)! onde temos que ter alguma informa¸c˜ao sobre a fun¸c˜ao que permita limitar sua derivada de ordem n + 1. ˜ CAP´ITULO 5. INTERPOLAC¸ AO POLINOMIAL 67 5.4 Escolha dos Pontos Uma das caracter´isticas da interpola¸c˜ao ´e que esta pode fornecer uma aproxima¸c˜ao local, sem a necessidade de usar todos os dados dispon´iveis. Como exemplo, vamos considerar a tabela abaixo x 0.2 0.34 0.4 0.52 0.6 0.72 f(x) 0.16 0.22 0.27 0.29 0.32 0.37 . Vamos achar uma aproxima¸c˜ao para f(0:44), usando um polinˆomio de grau 2. Neste caso, necessitamos de 3 pontos e o ideal ´e escolher aqueles que est˜ao mais pr´oximos do valor que desejamos aproximar. Logo a melhor escolha ser´a x0 =0:34;x1 =0:4e x2 =0:52. Isto se justifica pela f´ormula do erro, pois En(0:44)(0:44¡0:34)(0:44¡0:4)(0:44¡0:52)max jf(n+1)(x)| =0:00032 max jf(n+1)(x)j jj·j| x2[x0;x2] (n + 1)! x2[x0;x2] (n + 1)! Se tiv´essemos escolhido x0 =0:2e x2 =0:72, o erro estaria limitado por En(0:44)= 0:00268 max jf(n+1)(x)j: j| x2[x0;x2] (n + 1)! 5.5 Interpola¸c˜ao Inversa Considere o seguinte problema: Dada uma tabela de pontos (xk;fk) e um n´umero y . [f0;fn]. Desejamos achar o valor de x de tal forma que f(x)= y. Temos duas formas de resolver o problema: Obter o polinˆomio interpolador de f(x)e resolver a equa¸c˜ao pn(x)= y. Em geral a equa¸c˜ao pn(x)= y tem mais de uma solu¸c˜ao e se o grau do polinˆomio for maior que 2, n˜ao temos um procedimento anal´itico que determine as solu¸c˜oes. A outra forma de se achar x ´e fazer uma interpola¸c˜ao inversa. Se f(x) ´e invers´ivel num intervalo contendo y, ent˜ao interpolamos a fun¸c˜ao inversa, isto ´e consideramos o conjunto de dados x = f¡1(y) e achamos o polinˆomio interpolador de f¡1(y). A condi¸c˜ao para que a fun¸c˜ao seja invers´ivel ´e que esta seja mon´otona crescente ou decrescente. Em termos dos pontos tabelados isto significa que os pontos devem satisfazer f0 f1 > >fn ¢¢· ¢¢· Exemplo 5.5.1 Dada a tabela de pontos x 0.2 0.3 0.4 0.5 0.6 0.7 0.8 f(x) 0.587 0.809 0.951 1.000 0.951 0.809 0.587 ˜ CAP´ITULO 5. INTERPOLAC¸ AO POLINOMIAL 68 Vamos procurar uma aproxima¸c˜ao de x de tal forma que f(x)=0:9, usando uma interpola¸c˜ao quadr´atica na forma de Newton. Em primeiro lugar devemos determinar em que intervalo pode ocorrer f(x)=0:9. Neste exemplo temos duas possibilidades, para x . [0:3, 0:4] ou x . [0:6, 0:7]. Em segundo lugar devemos verificar se a fun¸c˜ao f(x) admite inversa. Para o primeiro caso temos que a fun¸c˜ao f(x) ´e crescente no intervalo [0:2, 0:5]. Logo esta admite inversa neste intervalo. No segundo caso a fun¸c˜ao admite inversa no intervalo [0:5, 0:8], pois esta ´e decrescente neste intervalo. Como desejamos uma interpola¸c˜ao quadr´atica temos que ter no m´inimo trˆes pontos e nos dois casos temos quatro pontos. Portanto ´e poss´ivel achar as duas aproxima¸c˜oes. Vamos nos concentrar no primeiro caso. Montando a tabela da fun¸c˜ao inversa temos y 0.587 0.809 0.951 1.000 f¡1(y) 0.2 0.3 0.4 0.5 Como desejamos um polinˆomio de grau 2, devemos escolher trˆes pontos e a melhor escolha s˜ao os pontos que est˜ao mais pr´oximos do valor a ser aproximado, y =0:9 (x0 =0:809;x1 =0:951 e x2 =1:000). Calculando as diferen¸cas divididas temos y f¡1 ordem 1 ordem 2 0.809 0.3 0.704 0.951 0.4 6.999 2.040 1.000 0.5 e conseq¨uentemente o polinˆomio ´e dado por p(y)=0:3+(y - 0:809)0:704 + (y - 0:951)(y - 0:809)6:999 =5:115 - 11:605y +6:994y 2 Portanto o valor de x tal que f(x)=0:9 ´e aproximado por x = f¡1(0:9) ˜ p(0:9) = 0:3315. 5.6 Observa¸c˜oes Finais Como no caso do ajuste de curvas, a interpola¸c˜ao aproxima uma fun¸c˜ao, sobre um conjunto de dados. Por´em a interpola¸c˜ao exige que os pontos sejam distintos. Fato que n˜ao precisa ocorrer com o ajuste de curvas. Al´em disso, o grau do polinˆomio interpolador depende da quantidade de pontos. Logo para um conjunto de 100 pontos teremos um polinˆomio de grau = 99, o que n˜ao ´e muito pr´atico se desejamos montar um modelo matem´atico. A interpola¸c˜ao ´e mais indicada para aproxima¸c˜oes quantitativas, enquanto que o ajuste de curvas ´e indicado para uma aproxima¸c˜oes qualitativas. Se desejamos saber a taxa de varia¸c˜ao de uma determinada fun¸c˜ao (ou seja a derivada), o mais indicado ´e o ajuste de curvas, pois na interpola¸c˜ao n˜ao temos garantia de que f0(x) ˜ p0n(x) (Veja figura 5.1). ˜ CAP´ITULO 5. INTERPOLAC¸ AO POLINOMIAL 69 Se a fun¸c˜ao que estamos interpolando, sobre n + 1 pontos ´e um polinˆomio de grau = n, ent˜ao o polinˆomio interpolador ´e a pr´opria f(x). Isto pode ser verificado pela f´ormula do erro, onde o termo f(n+1)(»)=0 8. . [x0;xn]. A interpola¸c˜ao ´e mais indicada para aproxima¸c˜oes quantitativas, enquanto que o ajuste de curvas ´e indicado para uma aproxima¸c˜ao qualitativa. A medida que aumentamos a quantidade de pontos num intervalo [a, b], ocorre o fenˆomeno de Runge, que ´e caracterizado em aumentar o erro nos pontos extremos do intervalo e melhorar a aproxima¸c˜ao nos pontos centrais. Isto ´e justificado pela f´ormula do erro, em que os pontos extremos do intervalo faz com que o fator (x - x0)(x - xn) seja grande. Desta ¢¢· forma, o polinˆomio interpolador n˜ao ´e indicado para extrapolar valores, isto ´e aproximar valores que n˜ao pertencem ao intervalo [x0;xn]. Abaixo apresentamos um exemplo implementado no MatLab, onde a figura 5.1 mostra o fenˆomeno de Runge. Figura 5.1: Fenˆomeno de Runge. ---pn(x), — f(x) % Diciplina de Calculo Numerico -Prof. J. E. Castilho % Forma de Lagrange do Pol. Interpolador % Interpola a funcao f(x)=1/(1+25x^2) nos pontos % x=[-1.0,-0.8,-0.6,...,0.6,0.8,1.0] % Calcula o polinomio e a funcao nos pontos ˜ CAP´ITULO 5. INTERPOLAC¸ AO POLINOMIAL 70 % xf=[-1.0,-0.98,-0.96,...,0.96,0.98,1.0] % Compara os graficos e mostra o fenomeno de Runge % clear; xf=-1:0.02:1; f=1./(1+25 *xf.^2); x=-1:0.2:1; fk=1./(1+25 *x.^2); % Forma de Lagrange n=size(xf); % pontos onde vamos calcular o polinomio m=size(x); % pontos de interpolacao for s=1:n(2) p(s)=0; for l=1:m(2) L=1; for k=1:l-1 L=L*(xf(s) -x(k))/(x(l)-x(k)); end; for k=l+1:m(2) L=L*(xf(s) -x(k))/(x(l)-x(k)); end; p(s)=p(s)+fk(l)*L; end; end; plot(xf,f,xf,p,':'); print -depsc fig51.eps 5.7 Exerc´icios Exerc´icio 5.1 A tabela abaixo fornece o n´umero de habitantes do Brasil (em milh˜oes) de 1900 a 1970. ano 1900 1920 1940 1950 1960 1970 Hab. 17.4 30.6 41.2 51.9 70.2 93.1 a) Usando o polinˆomio interpolador de grau 2, na forma de Lagrange, ache uma aproxima¸c˜ao para a popula¸c˜ao no ano de 1959. b) Usando interpola¸c˜ao quadr´atica na forma de Newton, estime, com o menor erro poss´ivel, em que ano a popula¸c˜ao ultrapassou os 50 milh˜oes. ˜ CAP´ITULO 5. INTERPOLAC¸ AO POLINOMIAL 71 c) Com os resultados obtidos no ´item a) podemos estimar a taxa de crescimento da popula¸c˜ao neste per´iodo? Justifique sua resposta. Exerc´icio 5.2 Considere a fun¸c˜ao f(x)= px e os pontos x0 = 16, x1 = 25 e x2 = 36. Com que precis˜ao podemos calcular p20, usando interpola¸c˜ao sobre estes pontos? Exerc´icio 5.3 Considere o problema de interpola¸c˜ao linear para f(x)= sen(x)+x, usando os pontos x0 e x1 = x0 + h. Mostre que jE1(x)j= h2=8. Exerc´icio 5.4 Dada a tabela abaixo. x -0.2 -0.1 0.1 0.15 0.35 f(x) 0.980 0.995 0.996 0.988 0.955 a) Quando poss´ivel, ache uma aproxima¸c˜ao para f(¡0:25) e f(0) , usando o polinˆomio interpolador na forma de Newton, com o menor erro poss´ivel. b) Se tiv´essemos usado o polinˆomio interpolador na forma de Lagrange sobre os mesmos pontos obter´iamos melhor resultado? Justifique. Exerc´icio 5.5 Num experimento de laborat´orio, uma rea¸c˜ao qu´imica liberou calor de acordo com as medi¸c˜oes mostradas na tabela abaixo hora 8:00 hr 9:00 hr 9:30 hr 10:00 hr Co 0 130 210 360 a) Determine uma fun¸c˜ao que dˆe a temperatura em fun¸c˜ao do tempo, de modo que os pontos tabelados sejam representados sem erro. b) Calcule a prov´avel temperatura ocorrida `as 9:45 hr. Cap´itulo 6 Integra¸c˜ao Num´erica -F´ormulas de Newton Cˆotes O objetivo deste cap´itulo ´e estudar esquemas num´ericos que aproxime a integral definida de uma fun¸c˜ao f(x) num intervalo [a, b]. A integra¸c˜ao num´erica ´e aplicada quando a primitiva da fun¸c˜ao n˜ao ´e conhecida ou quando s´o conhecemos a fun¸c˜ao f(x) num conjunto discreto de pontos. As f´ormulas de Newton-Cˆotes s˜ao baseadas na estrat´egia de aproximar a fun¸c˜ao f(x) por um polinˆomio interpolador e aproximamos a integral pela integral do polinˆomio. As aproxima¸c˜oes s˜ao do tipo n . b f(x)dx ˜ A0f(x0)+ A1f(x1)+ Anf(xn)= . Aif(xi) a ¢¢· i=0 onde os pontos s˜ao igualmente espa¸cados, isto ´e xk = x0 + kh, com h =(b - a)=n, eos coeficientes Ai s˜ao determinado pelo polinˆomio escolhido para aproximar a fun¸c˜ao f(x). 6.1 Regra do Trap´ezio A Regra do Trap´ezio ´e a f´ormula de Newton-Cˆotes que aproxima a fun¸c˜ao f(x) pelo polinˆomio interpolador de grau 1 sobre os pontos x0 = a e x1 = b. O polinˆomio interpolador de grau um, na forma de Newton ´e dado por p1(x)= f(x0)+ f[x0;x1](x - x0). Desta forma vamos aproximar a integral da fun¸c˜ao f(x) pela integral do polinˆomio, obtendo . b . x1 a f(x)dx ˜ x0 p1(x)dx . x1 = f(x0)+ f[x0;x1](x - x0)dx x0 72 CAP´ITULO 6. INTEGRAC¸ AO NUM ERICA -F ´ OTES 73 ˜´ ORMULAS DE NEWTON C ˆ 22 Ãx1 x0 . = f(x0)(x1 - x0)+ f[x0;x1] 2 - 2 - x0(x1 - x0) = f(x0)(x1 - x0)+ f(x1) - f(x0)(x1 - x0)2 x1 - x0 2 = (x1 - x0) (f(x0)+ f(x1)) 2 h =(f(x0)+ f(x1)); 2 onde h = A f´´ezio que tem f(x1)e f(x0) x1 - x0. ormula acima representa a area do trap´ como os valores das bases e h como o valor da altura. Na figura 6.1 temos uma representa¸c˜ao desta aproxima¸c˜ao. 6.2 C´alculo do Erro No cap´itulo de interpola¸c˜ao vimos que uma fun¸c˜ao f(x) pode ser representada por f(x)= pn(x)+ En(x), CAP´ITULO 6. ¸ AO NUM ERICA -F ´ ´ OTES 74 INTEGRAC˜ORMULAS DE NEWTON C ˆ onde pn(x) ´e o polinˆomio interpolador e En(x) o erro na interpola¸c˜ao definido no Teorema 5.3.1. Calculando a integral da fun¸c˜ao f(x) no intervalo [a, b], segue que . b . b . b f(x)dx = pn(x)dx + En(x)dx, aa a ou seja o erro na integra¸c˜ao ´e dado pela integra¸c˜ao do erro cometido na interpola¸c˜ao. No caso da Regra do Trap´ezio segue que o erro ´e dado por . b . b f00(») h3 ET = E1(x)dx =(x - x0)(x - x1) dx = - f00(»);. . [a, b] aa 2 12 Como o erro depende do ponto », que ´e desconhecido, na pr´atica usamos a estimativa h3 ET max f00(») jj= 12 »2[a;b] j| 2 Exemplo 6.2.1 Como exemplo, vamos considerar a fun¸c˜ao f(x)= e¡x, cuja a primitiva n˜ao tem uma forma anal´itica conhecida. Vamos aproximar a integral no intervalo [0, 1] usando a Regra do Trap´ezio. Desta forma temos que h =1 - 0=1 e segue que . 1 e¡x2 dx ˜ 1(e¡02 + e¡12 )=0:6839397 0 2 Para calcular o erro cometido temos que limitar a segunda derivada da fun¸c˜ao no intervalo [0, 1]. Sendo f00(x) = (4x 2 - 2)e¡x2 temos que nos extremos do intervalo vale jf00(0)| =2 e jf00(1)| =0:735759. Para calcular os pontos cr´iticos da f00(x), devemos derivar f00(x) e igualar a zero, obtendo, 2 s3 f000(x) = (12x - 8x 3)e¡x=0 x =0 ou x = § . 2 Como o ´unico ponto cr´itico pertencente ao intervalo ´e x =0 segue que max f00(»)=2 »2[a;b] j| Com isto temos que h3 1 ET max f00(»)= =0:166667 = 12 »2[a;b] j| 6 Note que esta estimativa do erro informa que a aproxima¸c˜ao obtida n˜ao tem garantida nem da primeira casa decimal como exata. Logo a solu¸c˜ao exata da integral est´a entre os valores 0:6839397 ± 0:166667. CAP´ITULO 6. INTEGRAC¸ AO NUM ERICA -F ´ OTES 75 ˜´ ORMULAS DE NEWTON C ˆ 6.3 Regra do Trap´ezio Repetida A Regra do Trap´ezio aproxima bem fun¸c˜oes suaves ( jf0(x)j. 1) e/ou em intervalos de integra¸c˜ao de amplitude pequena. Para intervalos de amplitude grande podemos fazer uma subdivis˜ao do intervalo de integra¸c˜ao [a, b] em n subintervalos de mesma amplitude e aplicamos a Regra do Trap´ezio em cada subintervalo (Ver Figura 6.1). Neste caso temos que Figura 6.1: Regra do Trap´ezio Repetida h = b - a e xk = x0 + kh com k =0, 1;:::;n n Aplicando a Regra do Trap´ezio em cada subintervalo [xk;xk+1] segue que . b hhh hh a f(x)dx ˜ 2(f0 + f1)+ 2(f1 + f2)+ 2(f2 + f3)+ ¢¢· + 2(fn¡2 + fn¡1)+ 2(fn¡1 + fn) h =[f0 + 2(f1 + f2 ++ fn¡1)+ fn] 2 ¢¢· O erro cometido na aproxima¸c˜ao ´e igual a soma dos erros cometidos em cada subintervalo, logo temos que o erro ´e da forma h3 max jETGj= n 12 »2[a;b] jf00(»)| CAP´ITULO 6. INTEGRAC¸ AO NUM ERICA -F ´ OTES 76 ˜´ ORMULAS DE NEWTON C ˆ 2 Exemplo 6.3.1 Considerando a fun¸c˜ao do exemplo anterior, f(x)= e¡xem [0, 1], vamos determinar o n´umero de subintervalos necess´arios para que a Regra do Trap´ezio Repetida forne¸ca uma aproxima¸c˜ao com pelo menos 3 casas decimais exatas. Para isto devemos ter que jETGj= 10¡4, logo h3 n max f00(»)= 10¡4 12 »2[0;1] j| Sendo h =(b - a)=n e que o m´aximo da segunda derivada da fun¸c˜ao em [0, 1] ´e 2 (Ver ex. anterior) segue que h3 11 n max f00(»)= 10¡4 n 2 = 10¡4 12 »2[0;1] j| . n3 12 1 2 . n= 6 × 10¡4 . n 2 = 1666:666 . n = 40:824 Devemos tomar no m´inimo n = 41 subintervalos para atingir a precis˜ao desejada. Abaixo apresentamos o uso da fun¸c˜ao trapz.m do MatLab que fornece a aproxima¸c˜ao da integral pela Regra do Trap´ezio. Usando 41 subintervalos temos a aproxima¸c˜ao It =0:746787657 % Diciplina de Calculo Numerico -Prof. J. E. Castilho % Exemplo do uso da funcao trapz(x,y) % Calcula a integral usando a regra do trapezio % para os pontos % x=[xo,x1,...,xn] % f=[fo,f1,...fn] % h=1/41; x=0:h:1; f=exp(-x.^2); It=trapz(x,f) 6.4 Regra de Simpson 1/3 Neste caso, usamos o polinˆomio de grau 2 que interpola a fun¸c˜ao f(x). Para isto necessitamos de trˆes pontos x0;x1 e x2. Como os pontos devem ser igualmente espa¸cados, tomamos h =(b - a)=2. Na figura 6.4 temos uma representa¸c˜ao desta aproxima¸c˜ao. Para obter a aproxima¸c˜ao da integral vamos considerar o polinˆomio interpolador na forma de Newton, p2(x)= f(x0)+(x - x0)f[x0;x1]+(x - x0)(x - x1)f[x0;x1;x2]. CAP´ITULO 6. INTEGRAC¸ AO NUM ERICA -F ´ OTES 77 ˜´ ORMULAS DE NEWTON C ˆ E com isto segue que a integral da fun¸c˜ao f(x) ´e aproximada por . b . x2 a f(x)dx ˜ x0 p2(x)dx h =(f(x0)+4f(x1)+ f(x2)); 3 De forma an´aloga a Regra do Trap´ezio, obtemos a f´ormula do erro, integrando o erro cometido na interpola¸c˜ao. Desta forma obtemos ES = . b E2(x)dx = . b (x - x0)(x - x1)(x - x2) f000(») dx = - h5 f(iv)(»);. . [a, b] aa 3! 90 Na pr´atica usamos a estimativa para o erro h5 ESmax f(iv)(»)jj= 90 »2[a;b] j| 2 Exemplo 6.4.1 Vamos considerar a fun¸c˜ao do exemplo anterior: f(x)= e¡xno intervalo [0, 1]. Para usar a Regra de Simpson temos que ter trˆes pontos. Desta forma tomamos b - a 1 - 01 h === 2 22 CAP´ITULO 6. INTEGRAC¸ AO NUM ERICA -F ´ OTES 78 ˜´ ORMULAS DE NEWTON C ˆ E com isto segue a aproxima¸c˜ao ´e dada por: . b 1 f(x)dx ˜ (f(0) + 4f(1=2) + f(1)) = 0:74718 a 6 A limita¸c˜ao do erro depende da limita¸c˜ao da quarta derivada da fun¸c˜ao no intervalo [0, 1], sendo: f(iv)(x) = (12 - 48x 2 + 16x 4)e¡x2 Calculando nos extremos do intervalo temos jf(iv)(0)| = 12 e jf(iv)(1)| =0:735759 Calculando os pontos cr´iticos temos f(v)(x)=(¡32x 5 + 160x 3 - 120x)e¡x2 =0 x =0 ou x = §0:958572 ou x = §2:02018 . Como o ponto x =0:958572 pertence ao intervalo [0,1] temos jf(iv)(0:958572)| =7:41948 Assim temos que max f(iv)(»)=12 »2[a;b] j| Obtemos desta forma que h5 ES max f(iv)(»)=0:00416667 == 90 »2[a;b] j| Desta forma podemos garantir que as duas primeiras casas decimais est˜ao corretas. 6.5 Regra de Simpson Repetida Podemos melhorar a aproxima¸c˜ao da mesma forma que fizemos com a Regra do Trap´ezio. Vamos dividir o intervalo de integra¸c˜ao em n subintervalos de mesma amplitude. Por´em devemos observar que a Regra de Simpson necessita de trˆes pontos, logo o subintervalo que aplicaremos a f´ormula ser´a da forma [xk;xk+2] o que implica que n deve ser um n´umero par. Neste caso temos que h = b - a e xk = x0 + kh k =0, 1;:::;n n Aplicando a Regra de Simpson em cada subintervalo [xk;xk+2] segue que . b hhh h a f(x)dx ˜ 3(f0 +4f1 + f2)+ 3(f2 +4f3 + f4)+ ¢¢· + 3(fn¡4 +4fn¡3 + fn¡2)+ 3(fn¡2 +4fn¡1 + f h =[f0 + 4(f1 + f3 ++ fn¡1)2(f2 + f4 ++ fn¡2)+ fn] 3 ¢¢· ¢¢· CAP´ITULO 6. INTEGRAC¸ AO NUM ERICA -F ´ OTES 79 ˜´ ORMULAS DE NEWTON C ˆ O erro cometido nesta aproxima¸c˜ao ´e a soma dos erros em cada subintervalo [xk;xk+2]e como temos n=2 subintervalos segue que: h5 = n 180 »2[a;b] jESR| max jf(iv)(»)| 2 Exemplo 6.5.1 Considerando a integral da fun¸c˜ao f(x)= e¡x, no intervalo [0, 1‘], vamos determinar o n´umero de subintervalos necess´arios para obtermos uma aproxima¸c˜ao com trˆes casas decimais exatas. Para isto limitamos o erro por ESR= 10¡4 . Sendo h =(b - a)=n e max»2[0;1] jf(iv)(»)| = 12 segue que j| h5 1 12 n max f(iv)(»)= 10¡4 n 5 180 »2[a;b] j| . n180 = 10¡4 1 . n4 = 15 × 10¡4 . n 4 = 666:66666 . n = 5:0813274 O menor valor de n que garante a precis˜ao ´e n =6. Note que a Regra de Simpson necessita de bem menos subintervalos que a Regra do Trap´ezio (n = 41). 6.6 Observa¸c˜oes Finais As f´ormulas de Newton-Cˆotes s˜ao obtidas aproximando a fun¸c˜ao por um polinˆomio interpolador. Quanto maior for o grau do polinˆomio, mais preciso ser´a o resultado obtido. Por´em a estrat´egia das f´ormulas repetidas permitem obter uma aproxima¸c˜ao com uma certa precis˜ao desejada, usando f´ormulas que s˜ao obtidas por polinˆomios de baixo grau. Pelas f´ormulas de erro podemos observar que a Regra do Trap´ezio ´e exata para polinˆomios de grau um, enquanto que a Regra de Simpson ´e exata para polinˆomios de grau trˆes. 6.7 Exerc´icios Exerc´icio 6.1 Calcule as integrais pela Regra do Trap´ezio e pela Regra de Simpson usando seis subintervalos. . 4 pxdx e . 0:6 dx 10 1+ x Exerc´icio 6.2 Calcule o valor de p com trˆes casas decimais exatas usando a rela¸c˜ao p . 1 dx = 4 0 1+ x2 CAP´ITULO 6. INTEGRAC¸ AO NUM ERICA -F ´ OTES 80 ˜´ ORMULAS DE NEWTON C ˆ Exerc´icio 6.3 Mostre que se f0(x) ´e cont´inua em [a, b] e que f00(x) > 0, 8x . [a, b], ent˜ao a aproxima¸c˜ao obtida pela Regra do Trap´ezio ´e maior que o valor exato da integral. Exerc´icio 6.4 Dada a tabela abaixo, calcule a integral :150:30f(x)dx com o menor erro 0 x 0.15 0.22 0.26 0.30 poss´ivel. f(x) 1.897 1.514 1.347 1.204 Exerc´icio 6.5 Baseado na Regra de Simpson, determine uma regra de integra¸c˜ao para a integral dupla . b . d f(x, y)dx dy ac Aplique a regra para calcular uma aproxima¸c˜ao para . 1 . 1 x 2 + y 2dx dy 00 Cap´itulo 7 Equa¸c˜oes Diferenciais Ordin´arias (P.V.I.) Muitos dos modelos matem´aticos nas ´areas de mecˆanica dos fluidos, fluxo de calor, vibra¸c˜oes, rea¸c˜oes qu´imicas s˜ao representados por equa¸c˜oes diferenciais. Em muitos casos, a teoria garante a existˆencia e unicidade da solu¸c˜ao, por´em nem sempre podemos obter a forma anal´itica desta solu¸c˜ao. Neste cap´itulo vamos nos concentrar em analisar esquemas num´ericos para solu¸c˜ao de Problemas de Valor Inicial (P.V.I.), para equa¸c˜oes diferenciais de primeira ordem. Isto ´e , achar a fun¸c˜ao y(x) tal que . y. = f(x, y) y(x0)= y0 A aproxima¸c˜ao de y(x) ´e calculada nos pontos x1;x2;x3;:::, onde xk = x0 + kh. O valor da fun¸c˜ao no ponto xk ´e aproximado por yk, que ´e obtido em fun¸c˜ao dos valores anteriores yk¡1;yk¡2;:::;y0. Com isto podemos classificar os m´etodos em duas classes: M´etodos de Passo Simples: S˜ao aqueles em que o c´alculo de yk depende apenas de yk¡1. M´etodos de Passo M´ultiplo: S˜ao aqueles em que o c´alculo de yk depende m-valores anteriores, yk¡1;yk¡2;:::;yk¡m. Neste caso dizemos que o m´etodo ´e de m-passos. 7.1 M´etodo Euler Dado o problema de valor inicial . y. = f(x, y) y(x0)= y0 Uma forma de aproximar a derivada de uma fun¸c˜ao no ponto x1 ´e dado por y0(x0) = lim y(x0 + h) - y(x0) y(x0 + h) - y(x0) (7.1) h!0 h ˜ h 81 ˜´ CAP´ITULO 7. EQUAC¸ OES DIFERENCIAIS ORDIN ARIAS (P.V.I.) Como x1 = x0 + h e pelo P.V.I. segue que y1 ¡h y0 = f(x0;y0) . y1 = y0 + hf(x0;y0) Com isto relacionamos o ponto y1 com y0, um valor dado no P.V.I. Assim obtemos uma aproxima¸c˜ao y(x1) ˜ y1. De forma an´aloga podemos obter y2 em fun¸c˜ao de y1, sendo que de uma forma geral teremos yk+1 = yk + hf(xk;yk) Este m´etodo ´e conhecido como M´etodo de Euler. Exemplo 7.1.1 Consideremos o seguinte problema de valor inicial . y. = x - 2y y(0) = 1 Neste caso temos que x0 =0 e y0 =1. Vamos usar o M´etodo de Euler para obter uma aproxima¸c˜ao para y(0:5), usando h =0:1. Desta forma temos y1 = y0 + h(x0 - 2y0)=1+0:1(0 - 2 * 1) = 0:8 ˜ y(x1)= y(0:1) y2 = y1 + h(x1 - 2y1)=0:8+0:1(0:1 - 2 * 0:8) = 0:65 ˜ y(x2)= y(0:2) y3 = y2 + h(x2 - 2y2)=0:65 + 0:1(0:2 - 2 * 0:65) = 0:54 ˜ y(x3)= y(0:3) y4 = y3 + h(x3 - 2y3)=0:54 + 0:1(0:3 - 2 * 0:54) = 0:462 ˜ y(x4)= y(0:4) y5 = y4 + h(x4 - 2y4)=0:462 + 0:1(0:4 - 2 * 0:462) = 0:4096 ˜ y(x5)= y(0:5) ˜´ CAP´ITULO 7. EQUAC¸ OES DIFERENCIAIS ORDIN ARIAS (P.V.I.) A solu¸c˜ao anal´itica do P.V.I. ´e dada por y(x) = (5e¡2x +2x - 1)=4. No gr´afico abaixo comparamos a solu¸c˜ao exata com os valores calculados pelo M´etodo de Euler. Note que em cada valor calculado o erro aumenta. Isto se deve porque cometemos um ”erro local” na aproxima¸c˜ao da derivada por (7.1) e este erro vai se acumulando a cada valor calculado. Na pr´oxima se¸c˜ao daremos uma forma geral do erro local. 7.2 M´etodos da S´erie de Taylor Vamos considerar o problema de valor inicial . y. = f(x, y) y(x0)= y0 Aplicando a s´erie de Taylor para y(x) no ponto xk, temos y(x)= y(xk)+ y0(xk) (x¡xk)+ y00(xk) (x¡xk)2 ++ y(n)(xk) (x¡xk)n + y(n+1)(») (x¡xk)n+1 1! 2! ¢¢· n!(n + 1)! Calculando no ponto xk+1 e considerando que xk+1 - xk = h temos que y0(xk) y00(xk) h2 y(n)(xk) hn y(n+1)(») hn+1 y(xk+1)= y(xk)+ h + ++ + (7.2) 1! 2! ¢¢· n!(n + 1)! Como y0(xk)= f(xk;yk) podemos relacionar as derivadas de ordem superior com as derivadas da fun¸c˜ao f(x, y). Como exemplo consideremos d y00(xk)= f(xk;yk) dx ˜´ CAP´ITULO 7. EQUAC¸ OES DIFERENCIAIS ORDIN ARIAS (P.V.I.) = fx + fyy. = fx + fyf y000(xk)= fy(fx + fyf)+ f2fyy +2ffxy + fxx Desta forma podemos obter uma aproxima¸c˜ao para o c´alculo do P.V.I., substituindo as rela¸c˜oes do tipo acima na s´erie de Taylor. Os m´etodos podem ser classificados de acordo com o termo de maior ordem que usamos na s´erie de Taylor, sendo Defini¸c˜ao 7.2.1 Dizemos que um m´etodo para a solu¸c˜ao de P.V.I. ´e de ordem n se este coincide com a s´erie de Taylor at´e o n-´esimo termo. O erro local cometido por esta aproxima¸c˜ao ser´a da forma y(n+1)(») hn+1Eloc(xk+1)= (n + 1)! . . [xk;xk+1] Como exemplo temos que o M´etodo de Euler ´e um m´etodo de 1o ordem, pois este ¯ coincide com a s´erie de Taylor at´e o primeiro termo, logo o erro local ´e dado por Eloc(xk+1)= y0. (») h2 . . [xk;xk+1] 2! Em geral, podemos determinar a ordem de um m´etodo pela f´ormula do erro. Se o erro depende da n-´esima derivada dizemos que o m´etodo ´e de ordem n - 1. Exemplo 7.2.1 Vamos utilizar o m´etodo de 2¯ o ordem para aproximar y(0:5), usando h = 0:1, para o P.V.I. . y. = x - 2y y(0) = 1 Neste caso temos que f(x, y)= x - 2y, logo segue que y0. = fx + fyf =1 - 2(x - 2y) Substituindo em (7.2) temos que y y(xk+1)= y(xk)+ 0(xk) h + y00(xk) h2 1! 2! h2 = y(xk)+ h(xk - 2y(xk))+ (1 - 2xk +4y(xk)) 2 Portanto o m´etodo ´e dado por h2 yk+1 = yk + h(xk - 2yk)+ (1 - 2xk +4yk) 2 ˜´ CAP´ITULO 7. EQUAC¸ OES DIFERENCIAIS ORDIN ARIAS (P.V.I.) Sendo x0 =0 e y0 =1 obtemos que h2 y1 = y0 + h(x0 - 2y0)+ (1 - 2x0 +4y0) 2 0:12 = 1+0:1(0 - 2 * 1)+ (1 - 2 * 0+4 * 1) = 0:825 2 h2 y2 = y1 + h(x1 - 2y1)+ (1 - 2x1 +4y1) 2 0:12 =0:825 + 0:1(0:1 - 2 * 0:825) + (1 - 2 * 0:1+4 * 0:825) = 0:6905 2 h2 y3 = y2 + h(x2 - 2y2)+ (1 - 2x2 +4y2) 2 0:12 =0:6905 + 0:1(0:2 - 2 * 0:6905) + (1 - 2 * 0:2+4 * 0:6905) = 0:58921 2 h2 y4 = y3 + h(x3 - 2y3)+ (1 - 2x3 +4y3) 2 0:12 =0:58921 + 0:1(0:3 - 2 * 0:58921) + (1 - 2 * 0:3+4 * 0:58921) = 0:515152 2 h2 y5 = y4 + h(x4 - 2y4)+ (1 - 2x4 +4y4) 2 0:12 =0:515152 + 0:1(0:4 - 2 * 0:515152) + (1 - 2 * 0:4+4 * 0:515152) = 0:463425 2 O gr´afico na figura 7.1 compara os resultados obtidos por este m´etodo, com os resultados obtidos pelo m´etodo de Euler. Os resultados obtidos pelo m´etodo de 2o ordem s˜ao mais precisos. Quanto maior a ¯ ordem do m´etodo melhor ser´a a aproxima¸c˜ao. A dificuldade em se tomar m´etodos de alta ordem ´e o c´alculo da rela¸c˜ao de y(n+1)(x)=[f(x, y)](n). 7.3 M´etodos de Runge-Kutta A estrat´egia dos m´etodos de Runge-Kutta ´e aproveitar as qualidades dos m´etodos da S´erie de Taylor (escolher a precis˜ao) sem ter que calcular as derivadas totais de f(x, y). O m´etodo de Runge-Kutta de 1¯ o ordem ´e o M´etodo de Euler, que coincide com o m´etodo da S´erie de Taylor de 1o ordem. Vamos determinar um m´etodo de 2o ordem, conhecido ¯¯ como m´etodo de Euler Melhorado ou m´etodo de Heun. Como o pr´oprio nome diz, a id´eia ´e modificar o m´etodo de Euler de tal forma que podemos melhorar a precis˜ao. Na figura 7.2-a temos a aproxima¸c˜ao yen+1 = yn + hf(xn;yn) que ´e obtida pela aproxima¸c˜ao do m´etodo de Euler. Montamos a reta L1 que tem coeficiente angular dado por y0(xn)= f(xn;yn). L1(x)= +(x - xn)y0= yn +(x - xn)f(xn;yn) ynn L1(xn+1)= yn +(xn+1 - xn)yn0 = yn + hf(xn;yn)= yen+1 ˜´ CAP´ITULO 7. EQUAC¸ OES DIFERENCIAIS ORDIN ARIAS (P.V.I.) Figura 7.1: M´etodo de Euler ’o’ —-M´etodo de 2¯ o ordem ’* ’ Montamos a reta L2 com coeficiente angular dado por f(xn+1, yen+1)= f(xn;yn + hy0)e passa pelo ponto P ( Ver figura 7.2-b ). L2(x)= yen+1 +(x - xn+1)f(xn+1, yen+1) Montamos a reta L0 que passa por P e tem como coeficiente angular a m´edia dos coeficientes angular de L1 e L2 (Ver figura 7.2-c). Finalmente a reta que passa pelo ponto (xn;yn) e ´e paralela a reta L0 tem a forma 1 L(x)= yn +(x - xn)(f(xn;yn)+ f(xn+1;yn + hy0)) 2n Calculando no ponto xn+1 temos 1 yn+1 = yn + h (f(xn;yn)+ f(xn+1;yn + hyn0 )) 2 Podemos observar que o valor de yn+1 (Ver figura 7.2-d) est´a mais pr´oximo do valor exato que o valor de yen+1. Este esquema num´erico ´e chamado de m´etodo de Euler Melhorado, onde uma estimativa do erro local ´e dado por h3 Eloc(xn)max y000(»)jj= 6 »2[xn;xn+1] j| ˜´ CAP´ITULO 7. EQUAC¸ OES DIFERENCIAIS ORDIN ARIAS (P.V.I.) Figura 7.2: M´etodo de Euler Melhorado Determinamos o m´etodo de Euler Melhorado por uma constru¸c˜ao geom´etrica. Tamb´em podemos obter uma demonstra¸c˜ao anal´itica. Desenvolvemos a s´erie de Taylor da fun¸c˜ao f(x, y) e calculando no ponto (xn+1;yn + hyn0 ). A express˜ao encontrada deve concordar com a S´erie de Taylor at´e a segunda ordem. Em geral um m´etodo de Runge-Kutta de segunda ordem ´e dado por yn+1 = yn + a1hf(xn;yn)+ a2hf(xn + b1h, yn + b2hyn0 ), ˜´ CAP´ITULO 7. EQUAC¸ OES DIFERENCIAIS ORDIN ARIAS (P.V.I.) onde 8 >a1 + a2 =1 . a2b1 =1=2 (7.3) >a2b2 =1=2 . O m´etodo de Euler Melhorado ´e obtido com a1 = a2 =1=2e b1 = b2 = 1. M´etodos de ordem superior s˜ao obtidos seguindo o mesmo procedimento. Abaixo apresentamos um m´etodo de 3o ¯ e 4o ¯ R-K 3o ¯ ordem: 214 yn+1 = yn + K1 + K2 + K3 939 K1 = hf(xn;yn) K2 = hf(xn + h=2;yn + K1=2) K3 = hf(xn +3h=4;yn +3K2=4) R-K 4o ¯ ordem: 1 yn+1 = yn +(K1 +2K2 +2K3 + K4) 6 K1 = hf(xn;yn) K2 = hf(xn + h=2;yn + K1=2) K3 = hf(xn + h=2;yn + K2=2) K4 = hf(xn + h, yn + K3) 7.4 M´etodos de Adans-Bashforth S˜ao m´etodos de passo m´ultiplos baseados na integra¸c˜ao num´erica. A estrat´egia ´e integrar a equa¸c˜ao diferencial no intervalo [xn;xn+1], isto ´e . xn+1 . xn+1 y0(x)dx = f(x, y(x))dx (7.4) xn xn ) . xn+1 y(xn+1)= y(xn)+ f(x, y(x))dx (7.5) xn Integra¸c˜ao Num´erica A integral sobre a fun¸c˜ao f(x, y) ´e aproximada pela integral de um polinˆomio interpolador que pode utilizar pontos que n˜ao pertencem ao intervalo de integra¸c˜ao. Dependendo da escolha dos pontos onde vamos aproximar a fun¸c˜ao f(x, y) os esquemas podem ser classificados como: Expl´icito: S˜ao obtidos quando utilizamos os pontos xn;xn¡1;:::;xn¡m para interpolar a fun¸c˜ao f(x, y). Impl´icito: S˜ao obtidos quando no conjunto de pontos, sobres os quais interpolamos a fun¸c˜ao f(x, y), temos o ponto xn+1. ˜´ CAP´ITULO 7. EQUAC¸ OES DIFERENCIAIS ORDIN ARIAS (P.V.I.) 7.4.1 M´etodos Expl´icitos Vamos considerar o caso em que a fun¸c˜ao f(x, y) ´e interpolada sobre os pontos (xn;fn) e(xn¡1;fn¡1), onde fn = f(xn;yn). Considerando o polinˆomio interpolador na forma de Newton temos f(x, y) ˜ p(x)= fn¡1 + f[xn¡1;xn](x - xn¡1) Integrando sobre o intervalo [xn;xn+1] temos . xn+1 . xn+1 p(x)dx = fn¡1 + f[xn¡1;xn](x - xn¡1)dx xn xn 22 Ãxn+1 xn . = fn¡1(xn+1 - xn)+ f[xn¡1;xn] - xn¡1xn+1 - + xn¡1xn 22 = hfn¡1 + fn - fn¡1 1 ³xn+12 - 2(xn - h)xn+1 - xn 2 + 2(xn - h)xn . h 2 = hfn¡1 + fn - fn¡1 1 ³xn+12 - 2xnxn+1 + xn 2 +2h(xn+1 - xn). h 2 = hfn¡1 + fn - fn¡1 1 ³(xn+1 - xn)2 +2h2. h 2 h = (2fn¡1 +3fn) 2 Substituindo a aproxima¸c˜ao da integral em (7.5) obtemos o seguinte m´etodo h yn+1 = yn + (2fn¡1 +3fn) 2 Este m´etodo ´e um m´etodo expl´icito, de passo dois. Isto significa que yn+1 depende de yn e yn¡1. Logo necessitamos de dois valores para iniciar o m´etodo: y0 que ´e dado no P.V.I.; y1 que deve ser obtido por um m´etodo de passo simples. De uma forma geral, os m´etodos de k-passos necessitam de k-valores iniciais que s˜ao obtidos por m´etodos de passo simples, de ordem igual ou superior a ordem do m´etodo utilizado. Obtemos o m´etodo, aproximando a fun¸c˜ao pelo polinˆomio interpolador de grau um. Assim o erro local, cometido por esta aproxima¸c˜ao ser´a . xn+1 . xn+1 f00(», y(»)) E1(x)dx =(x - xn¡1)(x - xn) dx xn xn 2! 5 = h3 y000(») . . [xn;xn+1] 12 Com isto temos a seguinte estimativa para o erro local 5 Eloc(xn+1)= h3 max y000(») j| 12 »2[xn;xn+1] j| ˜´ CAP´ITULO 7. EQUAC¸ OES DIFERENCIAIS ORDIN ARIAS (P.V.I.) Exemplo 7.4.1 Considere o seguinte P.V.I. . y. = ¡2xy y(0:5) = 1 Vamos achar uma aproxima¸c˜ao para y(1:1) pelo M´etodo de Adams-Bashforth expl´icito , de passo dois, usando h =0:2. Do P.V.I segue que y0 =1 e x0 =0:5. Este ´e um m´etodo de passo dois e vamos ter que calcular y1 por um m´etodo de passo simples que tenha ordem igual ou superior ao do M´etodo de Adams-Bashforth. Como o erro local depende da 3o derivada de y(x) temos que ¯ este m´etodo ´e de 2o ordem. Assim vamos utilizar o m´etodo de Euler Melhorado, que ´e um ¯ m´etodo de 2¯ o ordem. 1 y1 = y0 + h (f(x0;y0)+ f(x0;y0 + hy00 )) 2 0:2 =1+ (¡2 * 0:5 * 1+(¡2) * 0:5 * (1 + 0:2 * (¡2 * 0:5 * 1))) = 0:82 ˜ y(0:7) 2 Tendo o valor de y0 e y1 podemos iniciar o M´etodo de Adams-Bashforth, sendo h y2 = y1 + (2f0 +3f1) 2 0:2 =0:82+ (2 * (¡2) * 0:5 * 1+3 * (¡2) * (0:7) * 0:82) = 0:2756 ˜ y(0:9) 2 h y3 = y2 + (2f1 +3f2) 2 0:2 =0:2756 + (2 * (¡2) * 0:7 * 0:82 + 3 * (¡2) * (0:9) * 0:2756) = ¡0:102824 ˜ y(1:1) 2 Se tiv´essemos aproximado a fun¸c˜ao f(x, y) por um polinˆomio de grau 3, sobre os pontos (xn;fn), (xn¡1;fn¡1), (xn¡2;fn¡2), (xn¡3;fn¡3) obter´iamos o m´etodo de passo 4 e ordem 4 dado por h 251 yn+1 = yn+ [55fn - 59fn¡1 + 37fn¡2 - 9fn¡3] Eloc(xn+1)= h5 y(v)(») com . . [xn;xn+1] 24 720 Neste caso necessitamos de quatro valores iniciais, y0;y1;y2 e y3, que deve ser calculados por um m´etodo de passo simples de ordem maior ou igual a quatro (Ex. Runge-Kutta de 4o ¯ ordem). 7.4.2 M´etodos Impl´icitos Neste caso o ponto (xn+1;fn+1) ´e um dos ponto, onde a fun¸c˜ao f(x, y) ser´a interpolada . Vamos considerar o caso em que a fun¸c˜ao f(x, y) ´e interpolada sobre os pontos (xn;fn)e (xn+1;fn+1). Considerando o polinˆomio interpolador na forma de Newton temos f(x, y) ˜ p(x)fn + f[xn;xn+1](x - xn) CAP´ITULO 7. ¸ OES DIFERENCIAIS ORDIN ARIAS (P.V.I.) EQUAC˜´ 91 Integrando sobre o intervalo [xn;xn+1] temos . xn+1 . xn+1 p(x)dx = fn + f[xn;xn+1](x - xn)dx xn xn Ãxn+12 xn 2 . = fn(xn+1 - xn)+ f[xn;xn+1] - xnxn+1 - + xnxn 2 2 = hfn + fn+1 - fn 1 ³xn+12 - 2xnxn+1 + xn 2. h 2 = hfn + fn+1 - fn 1(xn+1 - xn)2 h 2 h =(fn + fn+1) 2 Substituindo a aproxima¸c˜ao da integral em (7.5) obtemos o seguinte m´etodo h yn+1 = yn +(fn + fn+1) 2 O erro local, cometido por esta aproxima¸c˜ao ser´a . xn+1 . xn+1 f00(», y(»)) E1(x)dx =(x - xn+1)(x - xn) dx xn xn 2! 1 = ¡h3 y000(») . . [xn;xn+1] 12Note que o c´alculo de yn+1 depende de fn+1 = f(xn;yn+1). Em geral a f(x, y) n˜ao permite que isolemos yn+1. Desta forma temos que usar um esquema Preditor-Corretor. Por um m´etodo expl´icito encontramos uma primeira aproxima¸c˜ao para yn+1 (Preditor) e este valor ser´a corrigido atrav´es do m´etodo impl´icito (Corretor). Exemplo 7.4.2 Vamos considerar o seguinte problema: . y= ¡2xy - y . 2 y(0) = 1 Usando h =0:1 vamos achar uma aproxima¸c˜ao para y(0:2), usando o esquema P : yn+1 = yn + hf(xn;yn) h C : yn+1 = yn + (2fn + fn+1) 2 Sendo x0 =0 e y0 =1 temos P : y1 = y0 + hf(x0;y0)=1+0:1(¡2 * 0 * 1 - 12)=0:9 C : y1 = y0 + h (2f0 + f1)=1+ 0:1 ³¡2 * 0 * 1 - 12 - 2 * 0:1 * 0:9 - 0:92. =0:9005 ˜ y(0:1) 2 2 P : y2 = y1 + hf(x1;y1)=0:9005 + 0:1(¡2 * 0:1 * 0:9005 - (0:9005)2)=0:8013 C : y2 = y1 + h (f1 + f2)=0:9005+ 0:1 ³¡2 * 0:1 * 0:9005 - (0:9005)2 - 2 * 0:2 * 0:8013 - 0:80132. 22 =0:8018 ˜ y(0:2) ˜´ CAP´ITULO 7. EQUAC¸ OES DIFERENCIAIS ORDIN ARIAS (P.V.I.) 92 7.5 Equa¸c˜oes de Ordem Superior Uma equa¸c˜ao de ordem m pode ser facilmente transformadas em um sistema de m equa¸c˜oes de primeira ordem. Este sistema pode ser visto como uma equa¸c˜ao vetorial de primeira ordem e assim poderemos aplicar qualquer um dos m´etodos analisados nas se¸c˜oes anteriores. Como exemplo vamos considerar o problema de terceira ordem dado por 8.. .. y00. = x2 + y2 - y. - 2y00x y(0) = 1 y0(0) = 2 y00(0) = 3 Para transformar num sistema de primeira ordem devemos usar vari´aveis auxiliares. Fazemos y. = wy0. = w. e fazemos y0. = w. = zy00. = w0. = z0. Com isto a equa¸c˜ao acima pode )) ser representada por: 8. <. . y. = w w. = z z. = x2 + y2 - w - 2zx y(0) = 1 w(0) = 2 z(0) = 3 O sistema acima pode ser escrito na forma matricial, onde y0 0B . CA . . B. . B. . C. . CA + . B. . C. 01 0 0 0 y w. z. y 00 1 w = 2 ¡1 ¡2x zx Y . Y X A Desta forma temos a equa¸c˜ao vetorial . Y . = AY + X Y (0) = Y0 onde Y0 = (1, 2, 3)T . Vamos aplicar o m´etodo de Euler na equa¸c˜ao acima obtendo Yn+1 = Yn + h(AnYn + Xn) ou seja 37 1C 1C 1. 1. . B. 1. 1. 0B . 0B . . = . + h . + Tomando h =0:1, vamos calcular uma aproxima¸c˜ao para y(0:2). Calculando Y1 ˜ Y (0:1) temos que . = . + h . + . . . B. . B. . 6. . B. . C. . B. . 7. . C. 01 0 0 yn+1 yn yn 00 1 0 wn+1 wn wn 2 ¡1 ¡2xn zn+1 zn yn zn xn . B. . 6. . B. . C. . B. . C. 01 0 0 y1 y0 y0 00 1 0 w1 w0 w0 z1 z0 y0 ¡1 ¡2x0 z0 x0 2 ˜´ CAP´ITULO 7. EQUAC¸ OES DIFERENCIAIS ORDIN ARIAS (P.V.I.) . = . +0:1 . + . = 37 1C 1C 1C . B. . B. . B. . 6. . B. . C. . B. . C. . B. . C. 1 01 0 1 0 1:2 y1 2 00 1 2 0 2:3 w1 02 3 3 2:9 z1 1 ¡1 ¡2 * 0 Calculando Y2 ˜ Y (0:2) temos . B@ . = . + h z2 z1 y1 1C 1C . B. . B. . 6. . B. . CA . + z1 x1 1C . B. 37 . CA . . 01 0 0 y2 y1 y1 00 1 0 w2 w1 w1 2 ¡1 ¡2x1 37 . = . +0:1 . + . = 1C Note que neste caso n˜ao estamos achando apenas o valor aproximado de y(0:2), mas tamb´em 1C ¼. 1C 1C . B. . B. . B. . 6. . B. . C. . B. . C. . B. . C. 1:2 01 0 1:2 0 1:430 y2 2:3 00 1 2:3 0 2:590 w2 0:12 2:9 1:2 ¡1 ¡2 * 0:1 2:9 2:757 z2 de y0(0:2) e y00(0:2), sendo . B. . B@ . C. y(0:2) 1:430 0(0:2) 2:590 y y00(0:2) 2:757 7.6 Exerc´icios Exerc´icio 7.1 Dado o P.V.I. . y. = x - y y(1) = 2:2 a) Considerando h =0:2, ache uma aproxima¸c˜ao para y(2:6), usando um m´etodo de segunda ordem. b) Se tiv´essemos usado o m´etodo de Euler, com o mesmo passo h, o resultado obtido seria mais preciso? Justifique. Exerc´icio 7.2 O esquema num´erico abaixo representa um m´etodo para solu¸c˜ao de E.D.O. h 251 yn+1 = yn + [9fn+1 + 19fn - 5fn¡1 + fn¡2] com Eloc(xn+1)= h5 y(5)(») 24 720 Classifique o m´etodo de acordo com o n´umero de passos, se este ´e impl´icito ou expl´icito, a ordem do m´etodo, procurando sempre dar justificativas as suas respostas. Exerc´icio 7.3 Determine o m´etodo de Runge-Kutta, de 2o ordem, obtido com a1 =0;a2 =1 ¯ e b1 = b2 =1=2 (ver (7.3)). Aplique o m´etodo para o P.V.I. 8>: <> y. = ¡2py y y(1) = 0:25 ˜ ´ CAP´ITULO 7. EQUAC¸ OES DIFERENCIAIS ORDIN ARIAS (P.V.I.) Exerc´icio 7.4 Considere o P.V.I. . y. = y y(0) = 1 a) Mostre que quando aplicamos o m´etodo de Euler Melhorado ao problema temos . h2 !n+1 yn+1 = 1+ h + 2 b) Comparando com a solu¸c˜ao exata do problema, vocˆe esperaria que o erro tivesse sempre o mesmo sinal? Justifique. Exerc´icio 7.5 Considere o P.V.I. abaixo. . y. = x - 2y y(0) = 1 a) Ache as aproxima¸c˜oes, para y(x) no intervalo [0, 2], usando h =0:2 e o esquema num´erico baixo P : yn+1 = yn + hf(xn + h=2;yn + h=2y0) h C ; yn+1 = yn +(fn+1 + fn) 2 b) Sabendo que a solu¸c˜ao exata ´e dada por y(x) = (5e¡2x +2x - 1)=4, plote (novo verbo?) um gr´afico com os valores obtidos pelo esquema Preditor e pelo esquema Corretor. Compare os resultados. Exerc´icio 7.6 Determine uma aproxima¸c˜ao para y(1) utilizando o m´etodo de Euler com h =0:1, para o P.V.I. abaixo. +2y =0 . y0. - 3y. y(0) = ¡1y0(0) = 0