Apostila Cáculo Numérico, Notas de estudo de Engenharia Informática

Baixe Apostila Cáculo Numérico e outras Notas de estudo em PDF para Engenharia Informática, somente na Docsity! PONTIFÍCIA UNIVERSIDADE CATÓLICA DO RIO GRANDE DO SUL FACULDADE DE MATEMÁTICA CÁLCULO NUMÉRICO Notas de Aula – Aplicações – Exercícios Eliete Biasotto Hauser afeifiz] EPI 8 Seul SESUIB/Q EJed UOWuaN / JeIUIOUIOA 0EHe]9d1B)UI + apeptur - gpuouMapueA ap zum Um ersuggo elouauoda SajpauM SELIGJEIG EIA, sagduna ap ajonhy EIMOUIOS BUENO FEUDÉEO TUE LSA0) SopeIpento SOUINIA S0p 0UB|I? Fglas-seneo) "opel 4” UOstus Igpoep-ssnea boLipunh opdeubejuI solzadei] + seda FEDJBA OuEMBjDA , soja SO)LGLINN SOpOJA -=»—— oq “hd SEMUS SEÍUBIANO sepessmy song saJtour sao Senha ap puIajsIS sepe quo Sega ap ogSeUS OquaLeMaIaIpuOY SeLBUIPO SISIDUBJBJIQ sagóenha | Cs Erabuny auesas vos 7 SOPRA hd aja SEMpOBIEA 3 SEMP "seda À | SejuapusasueaL à seougabIy sagáenha / Sezjpoo) belaassig jauaa STeIa e SIHISUBJALIQ Sagenba SOjuA SO: ego IO Op E3U UBE) OEMENAN - UG BP 58[404 > SOUIZ 50P BUGS] ORE ZIPALMON aQEN]NIA OJU0 3 ENMBJEIS E.B.Hauser – Cálculo Numérico Se a representação do real y em F não é exata, é necessário utilizar um arredondamento. Os tipos de arredondamento mais conhecidos são: - Arredondamento para número mais próximo de máquina (Oy). - Arredondamento por falta, ou truncamento (∇ y). - Arredondamento por excesso (∆ y) Em F, geralmente, as operações de adição e multiplicação não são comutativas, associativas e nem distributivas, pois numa série de operações aritméticas, o arredondamento é feito após cada operação. Ou seja, nem sempre as operações aritméticas válidas para os números reais são válidas em F. Isto influi na solução obtida através de um método numérico. Assim, métodos numéricos matematicamente equivalentes, podem fornecer resultados diferentes. Medidas de Exatidão Quando se aproxima um número real x por x*, o erro que resulta é x-x*. Define-se: erro absoluto: EA = *xx − e o erro relativo: ER = x xx *− para 0x ≠ . A fim de ver o tipo de situação que pode ocorrer um erro relativo de grande magnitude, vamos considerar a diferença entre os números a seguir, por exemplo: x = 0,3721478693 y = 0,3720230572 yx − = 0,0001248121=0,1248121 x 10-3 Se os cálculos forem feitos em F(10, 5, -499, 499) com arredondamento Ox: x*= 0,37215, y*= 0,37202 e x*-y*= 0,00013 = 0,13000 x 10-3 Assim, o erro relativo entre os dois resultados é grande: %4 3-10 x 0,1248121 3-10 x 13000,03-10 x 0,1248121 yx -y*)*x(yx ≈ − = − −− Na resolução de um problema o valor exato da solução x pode ser desconhecido. Podemos usar duas aproximações sucessivas de x, definindo: ⎥ ⎥ ⎦ ⎤ ⎢ ⎢ ⎣ ⎡ ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ + −+++−=+ 1 1log3,0)1,( kx kxkx kxkxDIGSE µ o qual expressa o número de dígitos significativos exatos de kx em relação a 1kx + . Aqui µ representa a unidade de arredondamento da máquina (µ = t1 2 1 −β se o arredondamen- to for Ox ). 5 1 - Teoria dos Erros Algoritmos Numéricos O Cálculo Numérico tem por objetivo estudar e aplicar algoritmos numéricos para a solução de problemas, visando o menor "custo" e confiabilidade do resultado (observar tempo de execução, número de operações aritméticas, quantidade de memória, propagação do erro, etc). Um bom algoritmo numérico deve se caracterizado por: a) Independência da máquina :a execução do programa pode ser realizada em diferentes máquinas. b) Inexistência de Erro Aritmético: problemas de overflow e underflow devem ser detecta- dos a priori c) Inexistência de Erro Lógico: previsão de tudo o que poderá ocorrer durante o processo (divisão por zero, por exemplo) d) Quantidade finita de cálculos. e) Existência de um critério de exatidão. f) O erro deve convergir para zero com precisão infinita. g) Eficiência: produz a resposta correta com economia Passos para a resolução de um problema: tipos de erros A resolução de problemas reais envolve diversas etapas: Estudaremos algoritmos numéricos a fim de obter uma solução numérica do problema, a qual, geralmente, difere da solução exata. Possíveis fontes de erro que geram essa diferença são: a) Simplificação do modelo matemático ( algumas variáveis envolvidas são desprezadas); b) Erro nos dados de entrada ( com frequência provindos de dados experimentais); c) Truncamento (por exemplo, substituição de um processo infinito por um finito e linearizações); d) Erro de arredondamento em aritmética de ponto flutuante. Curiosidades: Some disasters caused by numerical errors. http://ta.twi.tudelft.nl/users/vuik/wi211/disasters.html Aplicação O volume de uma esfera de raio r pode ser obtido de 3 3r4 V π = . Estimar o volume do hemisfério de raio “e” representado graficamente ao lado. Utilizar arredondamento para número mais próximo de máquina (0x) em F( 10 , 4 , -98 , 98). Problema Real Modelagem Matemática Solução 6 E.B.Hauser – Cálculo Numérico Exercícios 1) No sistema MapleV estimar 3.8e− utilizando ∑ ∞ = − =− 0i !i i)x(xe e ∑ ∞ = ==− 0i !i ix 1 xe 1xe com 26 termos cada e comparar com 3.8e− ≈0.2485168271x10-3. 2) Em F(10 , 3 , -98 , 98) e arredondamento por truncamento estimar p(2.73) se: a) 55.0x62x53x)x(p ++−= b) 55.0x)6x)5x(()x(p ++−= Em ambos os casos estimar o erro absoluto ao comparar com p(2.73) ≈0.11917x10-1. Obs: Estimar p(x) pelo algoritmo usual 01 2 2 3 3... 1 1)( axaxaxa nxna nxnaxp +++++ − −+= exige n adições e (n2+n)/2 multiplicações enquanto que o algoritmo de Horner { 0)1)2...)2)1 1 (((...()( axaxaxnaxnaxna n xp +++−+−+ − = requer n adições e n multiplicações. 3) Sejam A, B, C e D matrizes genéricas de ordem 10x20, 20x50, 50x1 e 1x100 respectivamente. Utilizando a propriedade associativa, pode-se determinar o produto matricial AxBxCxD de diversas formas. Qual das duas abaixo é mais eficiente? Porque? a) Ax(Bx(CxD)) b) (Ax(BxC))xD 4) Representar o número real x na base 2 usando 8 algarismos significativos? Essa representação é exata? a) x=0.6 b) x=13.25 c) x= 2.47 5) Determinar o cardinal , regiões de underflow e overflow e todos elementos reais de: a) F(2,3,-1,2) b) F(3,2,-1,2) c) F(2,2,-2,2) 6) Representar, se possível, os números abaixo em utilizando arredondamento por truncamento( x∇ )e arredondamento para número mais próximo de máquina (Ox) em F(10,5,-2,2). a) 3 b) 3 200 c) 3 2000 d) e e) 3000 1 f) 2 7) Considerando: ∑ = = 10 12 1 i i A a) Calcular o valor de A utilizando precisão infinita.. b) Utilizando arredondamento por truncamento ( x∇ ) em F(10,3,-98,98), estimar o valor de A somando da direita para esquerda e após somando da esquerda para a direita. Comparar os resultados. 7 E.B.Hauser – Cálculo Numérico 4) Se bia + é raiz de ( )xp de grau 2≥n , então ( )xp pode ser fatorado: ( ) ( ) ( )xqbaax2xxp 222 ⋅++−= onde o grau de ( )xq é 2−n . Ex.: a) ( ) 2x2xx2xxp 234 −++−= tem raízes i±± 1 ,1 . ( ) ( ) ( )1x2x2xxp 22 −⋅+−= Ex.: b) ( ) 10167 23 −+−= xxxxp tem raízes i±3 ,1 . ( ) ( ) ( )1x10x6xxp 2 −⋅+−= 5)Se ( )xp é de grau ímpar, então ( )xp possui ao menos uma raiz real. 6) Uma raiz 0x de ( )xp tem multiplicidade m se ( ) ( ) ( ) ( ) 0010"0'0 ===== − xpxpxpxp mK e ( ) 00 ≠xpm Ex.:1) 20 =x é raiz de multiplicidade 2 ( ) 8x6x2xp 23 +−= 2) 20 =x é raiz de multiplicidade 3 de ( ) 8465 234 −++−= xxxxxp 10 2 - Resolução de Equações Algébricas e Transcendentes 7) Valor numérico de um polinômio: para calcular, de forma usual, ( )ixp são necessárias n adições e ( ) 2 1+nn multiplicações. O Método de Horner faz esse cálculo com n adições e n multiplicações: ( ) 0121 parênteses 1 )))(((( axaxaxaxaxp nn n +++++= − − K 321 KK Ex.: ( ) 4x)3x)1x)2x)4x3((((4x3xx2x4x3xp 2345 +−−+=−+−−+= 8) Se ( )xp é de grau n , então existe único polinômio de grau ( )xqn ,1− , tal que ( ) ( ) ( ) ( )αα pxqxxp +⋅−= . Se α é raiz de ( )xp então ( ) 0=αp . É o algoritmo de Briot-Ruffini utilizado para Deflacionar Raízes. Ex. ( ) ( )( ) ( ) ( )( ) ( ) ( )( ) ( ) 1483pe14846x10x3x3 22pe26x5x2x2 01pe10x6x1x1 10x16x7xxp 2 3 2 2 2 1 23 −=−−+−+⇒−= =++−−⇒= =+−−⇒= −+−= α α α 2.1.1-Enumeração das Raízes Enumerar as raízes de ( )xp é dizer quantas são as raízes e se positivas, negativas ou complexas. Regra de Descartes ou Regras de Sinais “O número de raízes reais positivas de ( )xp é igual ao número de variações de sinais na seqüência dos coeficientes ou menor do que este por um número inteiro par, sendo uma raiz de multiplicidade “m” contada como “m” raízes e não sendo considerados os coeficientes nulos”. O número de raízes reais negativas é obtido aplicando a regra de Descartes a ( )xp − Regra de Huat Se ( ) 00p ≠ e para algum k, 11 2 +− ×≤ kkk aaa então ( )xp possui raízes complexas. Regra da Lacuna Se ( ) 00 ≠p e para algum K, 0 e 0 11 >×= +− kkk aaa , então ( )xp tem raízes complexas. 11 E.B.Hauser – Cálculo Numérico Ex.: ( ) 153 345 −−−+= xxxxxp • ( )xp pode ter: 1 raiz real positiva, 2 raízes negativas e 2 complexas; ou • 1 raiz real positiva, nenhuma negativa, e 4 complexas. 2.1.2-Localização das raízes de p(x) Localizar as raízes de ( )xp é determinar a região do plano que contém todas as raízes. Cota de Cauchy: Toda raiz α (real ou complexa) de ( )xp satisfaz βα ≤ . Onde 0,lim 0 == ∞→ xxkkβ e n k n n k n nn k n n k a a x a a x a a x a a x 0122111 ++++= −−−− + K Ex.: ( ) 13136,068,137,33 234 −+−+−= xxxxxp 4 1 23 1 )3136,068,137,33( +++=+ kkkk xxxx 00 =x e 4 9575796715,340048339019,3 9552764745,336165556530,2 59519059204,361069395791,2 09469665820,364735839615,1 69397307712,3557483314773,0 195 184 173 162 151 ≤∴ == == == == == α MM xx xx xx xx xx 12 2 - Resolução de Equações Algébricas e Transcendentes 1r 2r 3r 4r b) Cota de Cauchy: 0,11205,72 04 23 1 =+++=+ xxxx kkkk 46483007996,4 46481382609,481469532045,4 96478405782,467425618494,3 16472953921,460308011161,3 48211602868,1 17 164 153 142 1 = == == == = x xx xx xx x M M c) Separação de raízes ( ) 5,57617335495,3511515,8775,276 54321012345 −−−−− −−−−− xp x d) Cálculo da raiz ( )4,34r ∈ utilizando o método da bissecção. ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) M002874148,00352783203,313 20354003906,3;03515625,320152911151,020354003906,312 50356445312,3;03515625,30401045242,050356445312,311 0361328125,3;03515625,30897548754,00361328125,310 037109375,3;03515625,31891500395,0037109375,39 0390625,3;03515625,30095143543,003515625,38 039625,3;03125,33883184832,00390625,37 046875,3;03125,31900454168,1046875,36 0625,3;03125,3405335188,003125,35 0625,3;3871886353,20625,34 125,3;3660400391,9125,33 25,3;300390625,2525,32 35;39375,625,31 kIkxfkxk − − Obtemos 0354039062,312x4r =≈ com ( ) 096,413x,12xDIGSE ≅ 2.3.2-Método da Iteração Linear Consiste em escrever a equação ( ) 0=xf na forma ( )xGx = . Os pontos x* que satisfazem a condição ( )** xGx = são ditos pontos fixos de ( )xG e geometricamente representam os pontos de intersecção da reta xy = com a curva ( )xGy = . ( )xG é dita função de iteração do método. Inicia-se a iteração com um valor 0x próximo da raiz, e as outras aproximações são dadas por: ( ) K,2,1,0,1 ==+ ixGx ii 15 E.B.Hauser – Cálculo Numérico A seqüência de aproximação ix , converge para a solução x* da equação ( ) 0=xf sob certas condições . A construção de G não é única. A escolha de uma G apropriada é dita “problema do ponto fixo. Ex. 062 =−+ xx . ( ) 16, 1 6,6),6)6) 5432 2 1 −=+ =−−=−=−= x G x GxGcxGbxxGa Embora não seja preciso usar métodos numéricos para encontrar as duas raízes reais 2,3 21 =−= αα e da presente equação, nota-se que: i) Tomando 1G e 5,10 =x , a seqüência }{ ix não converge para 2. ( )ii xGx 11 =+ ( ) ( ) ( ) ( ) ( ) M 8342,12078822 46095276,3475 00390625,59)0625,8(6 0625,8)75,3(6 75,3)5,1(6 415 314 2 213 2 112 2 011 −== −== −=−−== −=−== =−== xGx xGx xGx xGx xGx ii) Tomando 2G e 5,10 =x , a seqüência }{ ix converge para 2. ( ) ( ) ( ) ( ) ( ) ( ) ( ) M 10000298018,2x 69998807918,1x 50004768183,2x 39980924992,1x 40076263645,2x 9694363804,1x 61213203435,25,16x 627 526 425 324 223 122 021 == == == == == == =−== xG xG xG xG xG xG xG Teorema da Convergência: Seja α uma raiz isolada de f em [ ]ba, . Se i) G e G’são contínuas em [ ]ba, ; ii) ( )b,ax,1)x('G ∈∀〈 ; iii) ( ) ,...2,1,0k,b,a)kx(G1kxe0x =∈=+∈Ι , então a seqüência { kx }, gerada por )(1 kxGkx =+ ,converge para α. 16 2 - Resolução de Equações Algébricas e Transcendentes Ex: Utilizando o método da iteração linear calcule a raiz de ( ) 3xexf x += , com ( ) 5, 1 ≥+kk xxDIGSE ( ) 33 33 00 x x xx eex exxexf −=−=⇒ −=⇒=+⇒= Seja ( ) ( ) ( ) ( ) 33,00' 24,01'; 3 1' 3 3 ≅ ≅−−= −= G GexG exG x x * G e G’ são continuas em [-1,0] e ( ) ]0,1[1' −∈∀< xxG . Logo , a seqüência gerada por 31 ix i ex −=+ converge para α ]0,1[−∈∀x . Seja 5,00 −=x M 772882595,0773204044,0 772884374,0771636903,0 772877469,0777723518,0 772904269,0754152577,0 772800243,0846481725,0 105 94 83 72 61 −=−= −=−= −=−= −=−= −=−= xx xx xx xx xx ( ) 34,5),( 000003188,0 109 9 ≅ = xxDIGSE exf *G não tem Maximo nem mínimo local em [0,1], testa-se então só os extremos. Características do Método da Iteração Linear: Não garante a convergência para toda função continua. Necessita do calculo de G’(x). Pode ocorrer dificuldade para encontrar G(x). A convergência é linear para raízes simples (a cada passo do método o erro é reduzido por um fato constante). A velocidade de convergência depende de ( )xG' , quanto menor este valor, maior será a convergência. 17 E.B.Hauser – Cálculo Numérico k kk k1k x 12 xlnx2xx + + −=+ 34 3 2 1 0 xx 426302751,0x 42699599,0x 42,0x 5,0x = = = = = Logo a aproximação para a raiz é 426302751,0=α com 9 dígitos significativos corretos. 2) Calcular a raiz ]4,3[4 ∈r do polinômio dado anteriormente: 11205,72)( 234 +−+= xxxxxp . 40)4(")4(0)3(")3( 0 =⇒>⋅<⋅ xppepp 201564 11205,72 23 234 1 −−+ −−−+ −=+ kkk kkkk kk xxx xxxx xx 56 5 4 3 2 1 0 03524990,3 03525211,3 03709653,3 08982331,3 36397059,3 4 xx x x x x x x = = = = = = = Obs: O Método de Newton pode divergir devido ao uso da tangente, oscilando indefinidamente. Uma aproximação para a raiz é 03524990,34 =r com 9 dígitos significativos e 5 iterações. 20 2 - Resolução de Equações Algébricas e Transcendentes Isto acontece quando: i) Não há raiz real ii) Ocorre simetria de ( )xf em torno de α iii) O valor inicial 0x está longe da raiz exata, fazendo que outra parte da função prenda a iteração. 2.3.4-Método da Secante Uma desvantagem do Método de Newton-Raphson é o calculo do valor numérico de ( )xf ' a cada iteração. O método da secante contorna este problema, substituindo a derivada pelo quociente das diferenças: 1 1 )()()(' − − − − ≅ kk kk k xx xfxf xf onde kx e 1−kx são duas aproximações para a raiz de ( )xf . A formula iterativa é: )()( )()( 1 1 1 − − + − ⋅− = kk kkk k xfxf xfxx x Geometricamente, a partir das aproximações para a raiz de kx e 1+kx , o ponto 1+kx é dado pela abscissa do ponto de intersecção do eixo Ox e da reta secante que passa por ))(,())(,( 11 kkkk xfxexfx −− . Características do método da secante: Garante a convergência para toda função continua Pode divergir se )()( 1−≅ kk xfxf Convergência mais lenta que o Método de Newton e mais rápida que Bissecção e Iteração Linear. São necessárias duas aproximações da raiz a cada iteração. Ex: )0,1(e21x172x53x)x(p −∈++−= α )21175()21175( )21175()( 1 2 1 3 1 23 23 1 1 ++−⋅++− ++−⋅− −= −−− − + kkkkkk kkkkk kk xxxxxx xxxxx xx 67 6 5 4 3 2 1 10 629321148566,0 89321148568,0 069321123947,0 869317876515,0 599601423487,0 898888888888,0 1,1 xx x x x x x x xx = −= −= −= −= −= −= =−= 21 E.B.Hauser – Cálculo Numérico Exercícios 1) Uma partícula é arremessada verticalmente, a partir do solo, com uma velocidade inicial ov .Desprezando a resistência do ar, supomos que a posição p partícula é dada por: 2 2 tgtov)t(d −= , onde g é aceleração da gravidade. Determinar a altura máxima atingida pela partícula e o instante em que ocorreu. 2) Uma corrente oscilante num circuito elétrico é descrita por I = )t2sen(te9 π− , t em segundos. Determinar todos os valores positivos de t para os quais I = 3.5. 3) A concentração de uma bactéria poluente num lago é descrita por C = t075.0e5.2t5.1e70 −+− Encontrar o tempo para que a concentração seja reduzida para nove. 4) O deslocamento de uma estrutura está definido pela seguinte equação D = wtcoskte8 − onde k = 0.5 e w = 3. a) Determinar graficamente uma estimativa inicial do tempo necessário para o deslocamento decrescer para 4. b) Usar o método de Newton-raphson para determinar essa raiz 5) Enumerar, localizar, separar e calcular (via Newton-Raphson e/ou Bissecção ), se possível, todas as raízes dos polinomios tendo como critério de parada DIGSE (xk , xk+1) ≥ 5. No caso de raízes múltiplas, determinar a multiplicidade da raiz e calcular as demais utilizando deflação. a) p(x) = x x x3 22 3 5− + − b) p(x) = x5 - 15,5x4 +77,5x3 - 155x2 +124x -31 c) p(x) = x x x x4 3 2121 2247 15043 34300− + − + d) p(x) = x x x x4 3 2115 1575 7625 12500− + − + e) p(x) = x x x x4 3 23 337 168 0 3136− + − +. . . f) p(x) = x4-11,101x3+11,1111x2-1,011x+0,001 g) p(x) = x9- 4x8 + 3,9x7 +3,02x6 - 5,5295x5 - 0,84732x4 + 2,83536x3 + 0,37152x2 -0,59616x - 0,15552 h) p(x) = x3 - 2081,93x2 +1424,64x- 3,19728 i) p(x) = x4 + 1,98x3 +1,1424x2 +0,162602x - 0,00174225 6) Localizar graficamente e calcular ( via Newton-Raphson e/ou Método da Iteração Linear) todas as raízes, com DIGSE(xk , xk+1) ≥ 5: 22 3. Sistemas de Equações Lineares O sistema com n equações lineares e n variáveis bxaxaxaxa bxaxaxaxa bxaxaxaxa nnnn33n22n11n 2nn2323222221 1nn1313212111 =++++ =++++ =++++ L MMMMM L L pode ser representado matricialmente por BAX = , onde ⎥ ⎥ ⎥ ⎥ ⎦ ⎤ ⎢ ⎢ ⎢ ⎢ ⎣ ⎡ = aaa aaa aaa A nn2n1n n22221 n11211 K MMM K K , ⎥ ⎥ ⎥ ⎥ ⎦ ⎤ ⎢ ⎢ ⎢ ⎢ ⎣ ⎡ = x x x X n 2 1 M , ⎥ ⎥ ⎥ ⎥ ⎦ ⎤ ⎢ ⎢ ⎢ ⎢ ⎣ ⎡ = b b b B n 2 1 M e A é a matriz dos coeficientes, X é o vetor das incógnitas e B o vetor dos termos independentes. 3.1- Introdução à problemática de sistemas Um SEL pode ser mal condicionado, bem condicionado ou sem solução. Um sistema é “mal condicionado” se uma pequena alteração nos dados pode provocar uma grande alteração na solução. Por exemplo: (a) ⎩ ⎨ ⎧ =+ =+ 5yx 95,4y98,0x tem solução exata ⎟⎟ ⎠ ⎞ ⎜⎜ ⎝ ⎛ = 5,2 5,2 x (b) ⎩ ⎨ ⎧ =+ =+ 5yx 95,4y99,0x tem solução exata ⎟⎟ ⎠ ⎞ ⎜⎜ ⎝ ⎛ = 0,5 0,0 x Uma alteração de 1% (0,01 no coeficiente 0,98) ocasionou uma alteração de 100% na solução. No caso de um sistema linear de ordem 2, cada equação representa uma reta. Resolver o sistema significa determinar a intersecção das duas retas. Três casos são possíveis: Bem condicionado Não tem solução. Mal condicionado 0det ≠ 0det = 0det ≅ retas concorrentes retas paralelas (perturbação em 2ℜ ) R1 R2 R2 R1 R2 R1 25 E.B.Hauser – Cálculo Numérico Exemplo2: (a) ⎩ ⎨ ⎧ =+ =+ 503,25y501,7x5,1 17y5x tem solução exata: ⎟⎟ ⎠ ⎞ ⎜⎜ ⎝ ⎛ = 3 2 x (b) ⎩ ⎨ ⎧ =+ =+ 501,25y501,7x5,1 17y5x tem solução: ⎟⎟ ⎠ ⎞ ⎜⎜ ⎝ ⎛ = 1 12 x 3.2-Medidas de Condicionamento O determinante normalizado da matriz dos coeficientes A é dado por 2 kn 2 2k 2 1kk n21 aaaondeAdetANORM L L ++== α ααα , para k = 1, 2, ..., n ( )1,1ANorm −∈ e quanto mais afastado de 1± (isto é, quanto mais próximo de zero) mais mal condicionado é a matriz A. Retomando o Ex2: 2) ⎟⎟ ⎠ ⎞ ⎜⎜ ⎝ ⎛ = 501,75,1 51 A ( ) ( ) 586495509853,7501,75,1 90990195135,5251 2 1 22 2 2 1 1 =+= =+= α α 38740000256377,0 0050000128,39 001,0001,0Anorm 001,05,15501,7Adet 21 == ⋅ = =⋅−= αα Agora, como pode ser medido o condicionamento de um sistema linear? Dado um SEL BAX = , o seu número de condicionamento é dado por: 1AA)A(Cond −= . Quanto maior for )A(Cond , mais sensível será o sistema linear. Utilizamos ∑== =≤≤ ∞ n 1i ij ni1 amaxAA , a norma do máximo das linhas. Ex: ⎟⎟ ⎠ ⎞ ⎜⎜ ⎝ ⎛ = 501,75,1 51 A , ⎟⎟ ⎠ ⎞ ⎜⎜ ⎝ ⎛ − − =− 10001500 50007501 A 1 51 1 101501,112521AA)A(Cond 12501A,001,9A ⋅≅== == − − 26 Sistemas de Equaçõe Lineares 3.3-Método de Resolução de Sistemas Métodos Diretos: A solução exata é obtida realizando-se um número finito de operações aritméticas em ℜ (com precisão infinita): Eliminação de Gauss e Gauss Jordan. Métodos Iterativos: A solução x é obtida como limite de uma seqüência de aproximações sucessivas x1, x2, ... . Método de Eliminação de Gauss Algoritmo básico de Gauss: A solução de BAX = é calculada em duas etapas: 1º- Triangularização : Mediante operações elementares nas linhas, a matriz A é transformada numa matriz triangular superior. Algoritmo: para 1n,,1k −= K (indica a linha do pivô) para n,,1ki K+= (indica a linha a transformar de A) kk ik a am −= 0aik = para nkj ,,1K+= (indica a coluna a transformar da linha i) amaa kjijij ⋅+= bmbb kii ⋅−= Obs.: Se 0aii = é necessário trocar de linha, se possível. 2º- Retro-substituição: A triangularização transforma o sistema BAX = , no sistema equivalente: dxc dxcxc dxcxcxc dxcxcxcxc nnnn 3nn3333 2nn2323322 1nn1313212111 = =++ =+++ =++++ LLLLLLL L L L cuja solução é dada por: ( ) ( ) 11 nn13132121 1 1n,1n nn,1n1n 1n nn n n a xcxcxcdx , a xcd x, c dx +++− = − == −− −− − L Teorema: O método de Gauss produz sempre a solução exata do sistema BAX = (utilizando precisão infinita) se 0det ≠A e as linhas de A forem permutadas quando 0aii = . 27 E.B.Hauser – Cálculo Numérico Exemplo: Resolver o sistema abaixo por Gauss- Jacobi e Gauss-Seidel. Critério de Parada: erro absoluto da solução calculada for menor que 10-3. ⎪ ⎪ ⎩ ⎪ ⎪ ⎨ ⎧ − − =+ =+ =+− =−+ 5,1 1 3 2 xx5,0 xx2 x2x xx4x 31 41 43 421 Reordenamos a fim de satisfazer ao critério de convergência. ⎪ ⎪ ⎩ ⎪ ⎪ ⎨ ⎧ − − =+− =+ =−+ =+ 3 5,1 2 1 x2x xx5,0 xx4x xx2 43 31 421 41 12 5,01 114 12 −> > −+> > Como a matriz dos coeficientes , após a reordenação, é diagonalmente dominante, então a aplicação dos métodos de Gauss-Jacobi e Gauss-Seidel produzirá uma seqüência ( )( )kX convergente para a solução exata. Gauss-Jacobi Fórmulas iterativas: ( ) ( ) ( ) ( ) ( ) ( ) ( )( ) ( ) ( )k 3 )k( 31k 4 k 1 1k 3 k 4 k 1 )k( 4 )k( 11k 2 k 4 )k( 41k 1 x5,05,1 2 )x3( x x5,05,1x x25,0x25,05,0 4 )xx2( x x5,05,0 2 x1 x +−= +− = −= +−−= +−− = −= ⎟ ⎠ ⎞ ⎜ ⎝ ⎛ − = + + + + Aproximação inicial: ( ) 0X 0 = . 30 Sistemas de Equaçõe Lineares k x1 x2 x3 x4 0 0 0 0 0 1 0,5 0,5 1,5 -1,5 2 1,25 -1,0 1,25 -0,75 3 0,875 -1,0 0,875 -0,875 4 0,9375 -0,9375 1,0625 -10625 5 1,03125 -1,0000 1,03125 -0,96875 6 0,984375 -1,0000 0,984375 -0,984375 7 0,9921875 -0,9921875 1,0078125 -1,0078125 8 1,00390625 -1,0000 1,00390625 0,9960375 9 0,998046875 -1,0000 0,998046875 -0,998046875 10 0,9990234375 -0,999023437 1,000976563 -1,000976563 11 1.000488282 -1.0000000 1.000488281 -0,999511718 : : : : : 40 1 -1 1 -1 na 12a. iteração consegue-se ( ) 51112 10xx −<− Gauss-Seidel Fórmulas iterativas: ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( )1k 3 1k 4 1k 1 1k 3 k 4 1k 1 1k 2 k 4 1k 1 x5,05,1x x5,05,1x x25,0x25,05,0x x5,05,0x ++ ++ ++ + +−= −= +−−= −= k x1 x2 x3 x4 0 0 0 0 0 1 0,5 -0,625 1,25 -0,875 2 0,9375 0,953125 1,03125 -0,984375 3 0,9921875 -0,994140625 1,00390625 0,99806875 4 0,9990234375 -0,9992675781 1,000488281 -0,999755895 5 0,999877929 -0,999908447 1,000061035 -0,999969482 : : : : : 12 1 -1 1 -1 ( ) ( ) 545 10xx −<− 31 E.B.Hauser – Cálculo Numérico EXERCÍCIOS 1) Uma consideração importante no estudo da transferência de calor é a de se determinar a distribuição de temperatura numa placa, quando a temperatura nas bordas é conhecida. Supomos que a placa da figura represente um seção transversal de uma barra de metal, com fluxo de calor desprezível na direção perpendicular à placa. Sejam T1, T2 , ..., T6 as temperaturas nos seis vértices interiores do reticulado da figura. A temperatura num vértice é aproximadamente igual à média dos quatro vértices vizinhos mais próximos (à esquerda, acima, à direita e abaixo). Por exemplo: T1 = ( 10 + 20 + T2 + T4 )/4 ou 4T1 = 10 + 20 + T2 + T4 a) Escrever o sistema de seis equações, cuja solução fornece estimativas para as temperaturas T1, T2 , ..., T6. b) Resolver o sistema utilizando o sistema MapleV. 2) Num experimento num túnel de vento, a força sobre um projétil devido à resistência do ar foi medida para velocidades diferentes. 119 10 3.74 8 6.39 6 8.14 4 90.2 2 0 0 força velocidade Expressar a força como função da velocidade aproximando-a a um polinômio de quinto grau: 5 5 4 4 3 3 2 21o vavavavavaa)v(f +++++= Verificar a validade de )(vf encontrada e obter uma estimativa para a força sobre o projétil quando ele estiver se deslocando a uma velocidade de 1, 3, 5, 7 e 9 unidades de velocidade. 3) Considere o sistema (#) ⎪ ⎩ ⎪ ⎨ ⎧ =++ =++ =++ 78,0x2,0x25,0x33,0 08,1x25,0x33,0x5,0 83,1x33,0x5,0x 321 321 321 . a) (#) é bem ou mal condicionado? Porque? O que isso significa? b) Resolver (#) pelo método de Gauss sem pivotamento. 1 20o 20o 20o 20o 20o 20o 40o 40o 10o 10o 1 2 3 4 5 6 32 Sistemas de Equaçõe Lineares validade := [0, 2.899999998, 14.80000000, 39.60000000, 74.29999994, 119.0000000] > estimativas:=[f(1),f(3),f(5),f(7),f(9)]; estimativas := [1.111718750, 7.202343750, 25.73046873, 55.89609367, 94.9992188] 3-a) Mal condicionado. NORM A ≅ 0,000181. Uma pequena perturbação nos dados de entrada pode causar uma grande alteração na solução. b) A solução exata é X=[1 1 1]T 4 - NORM A ≅ 0,00257, a solução exata é X=[1 1 1 1 1 1 1]T 5- a)X=[-21,86188 11,46568 2,376447 -8,514801 0,7475478 -15,50981 18,08498]T b) solução exata X=[1 0 2]T c) solução exata X=[1 2 3 -1]T 6-i) a) detA=0 e o sistema e incompatível b) detA = 0 e o sistema tem infinitas soluções dadas por x=z e y=-z ii) A não é matriz Diagonal Dominante e Gauss-Jacobi e Gauss-Seidel a) converge (oscilando) para a solução exata [1.5 1.5]T b) diverge 7 - a) X(10) = [-0.930569 1.901519 1.359500 -1.906078]T X(35) = [-1.076888 1.990028 1.474477 -1.906078]T b) solução exata X = [2 1 4]T 8 - A solução exata é x =1 y= 2 z = 3. 9 - soluções exatas: a) X = [3 -4 2]T b) X = [2 -7 4]T c) X = [-3 2 1]T 10 - a) detA = -55 b) det A = 706.875 11 - R = [0.12 0.24 0.25]T 35 4.Interpolação Polinomial 4.1 Objetivo: Dado um conjunto de n+1 pontos distintos ))x(f,x( ii , n,...,1,0i = , queremos determinar o polinômio p(x) talque ),x(f)x(p ii = e para os demais pontos do intervalo )x(f)x(p ≅ , [ ]no x;xx∈∀ . O polinômio p(x) é dito polinômio aproximante ou interpolador de f(x) no intervalo [ ]nxox , . 4.2 Aplicações • Obter uma expressão analítica aproximada de uma função que é conhecida em apenas um número finito de pontos; • Avaliar a função num ponto não tabelado [ ]nxoxx ;*∈ ; • Determinar aproximações para ∫ n x ox dxxf )( , substituindo a f(x) pelo polinômio interpolador; • Calcular uma aproximação para )(' xf para [ ]nxoxx ;∈ , substituindo f(x) por p(x). 4.2 Existência e Unicidade da Solução Dados ℜ∈ix e ,)x(f i ℜ∈ n,...,1,0i = , procuramos nP)x(p ∈ tal que iii x),x(f)x(p ∀= . Seja ∑ = =++++= n k kxka nxnaxaxaoaxp 0 ...221)( . 36 E.B.Hauser – Cálculo Numérico Então de =)x(p i ∑ = = n 0k i k ik ),x(fxa obtemos: ⎪ ⎪ ⎪ ⎩ ⎪⎪ ⎪ ⎨ ⎧ =++++ =++++ =++++ n n nn 2 n2n10 1 n 1n 2 12110 0 n 0n 2 02010 fxa...xaxaa .............................................. fxa...xaxaa fxa...xaxaa , o qual representa um sistema linear de ordem n+1 onde as n+1 incógnitas são os ,ka n,0k = e a matriz dos coeficientes é dada por: A = ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ ⎤ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ ⎡ n n 2 nn n 2 2 22 n 1 2 11 n 0 2 00 x..xx1 ...... ...... ...... x..xx1 x..xx1 x..xx1 De acordo com Vandermonde, det A = ⎥ ⎥ ⎦ ⎤ ⎢ ⎢ ⎣ ⎡ −∏∏ += − = n 1ij ij 1n 0i )xx( . Como os pontos são distintos, a diferença ixjx − será sempre diferente de zero, e então detA ≠ 0. Portanto o polinômio interpolador existe e também é único. 4.3Polinômio Interpolador de Newton Para Diferenças Finitas Ascendentes Dados ( )ii y,x , )x(fy ii = i=0,1,2,...,n satisfazendo hxx....xxxx 1nn1201 =−==−=− − , isto é , hxx i1i =−+ . Para k= 1, 2, ...n, e i= 0, 1, 2, ..., n-k , a diferença finita de ordem k é dada por i 1k 1i 1k i k yyy −+ − −= ∆∆∆ . E, o polinômio interpolador de Newton para diferenças finitas ascendentes é dado por : o n n 1n1o o 2 2 1o 0 o o y h!n )xx()xx)(xx( y h!2 )xx)(xx( y h )xx( y)x(p ∆∆∆ − −−− ++ −− + − += L L 37 4-Interpolação Polinomial Exemplo2: A velocidade do som na água varia com a temperatura conforme tabela: 532,1 110 538,1 4,104 544,1 9,98 548,1 3,93 552,1 86 )s/m(velocidade ) C(atemperatur o Solução: 1) Cálculo das diferenças divididas x y 0y∆ 0 2y∆ 0 3y∆ 0 4y∆ 86 1.552 -0.00054794 -0.00001289 -0.114 x 10-5 0.136 x 10-6 93.3 1.548 -0.00071428 -0.00003393 0.213 x 10-5 --------------- 98.9 1.544 -0.00109090 0.175 x 10-5 --------------- --------------- 104.4 1.538 -0.00107142 --------------- --------------- --------------- 110 1.532 --------------- --------------- --------------- --------------- 2) Construção do polinômio: )9.98x)(3.93x)(86x()00001289.0)(3.93x)(86x()00054794.0)(86x(552.1)x(p −−−+−−−+−−+= (-0.114 x 10 5− )+ )4.104x)(9.98x)(3.93x)(86x( −−−− (0.136 x 10 6− ) Simplificando a expressão , encontramos 5036369194.0x830077947032.0x99260000534327.0x10840.13666902p(x) 234-6 −+−×= 3) Verificação da validade de p(x) calculado no item 2: 1.55199999)86(p = ; 548.1)3.93(p = ; 544.1)9.98(p = ; 537999999.1)4.104(p = ; 532.1)110(p = 4) Estimativa da velocidade do som quando a temperatura da água atinge 100 C0 , é 1.54293925)100(p ≅ m/s 6) Análise gráfica do problema: 40 E.B.Hauser – Cálculo Numérico 4.5 Erro de Truncamento [ ]nxxnfn xxE ,0),( )1( )!1( )()( ∈+ + = ξξϕ com i 1n n101n 0 1n n10 y)*xx)...(xx)(xx( h)!1n( f)xx)...(xx)(xx()x( + + + −−−= + −−−= ∆ ∆ ϕ , pois i k ki k f h!.k 1y ∆∆ = , para h constante. 4.6 Aplicações utilizando o sistema Maple APLICAÇAO 1 Cinquenta animais de uma espécie ameaçada de extinção foram colocados numa reserva e um controle populacional mostrou os dados: 69 13 73 10 76 7 77 5 73 3 60 1 50 0 animaisdequantidade )anos(t a) Determinar a matriz de Vandermonde para o problema e determinar o valor do respectivo Número de Condicionamento (Cond e Determinante Normalizado). O que podemos concluir? b) Determinar o polinômio interpolador utilizando todos os dados tabelados. c) Verificar a validade do modelo encontrado. d) Estimar a a população no quarto ano. É possível estimar a população no décimo quinto ano utilizando o polinômio determinado no ítem b. e) Em que ano essa população animal atingiu seu máximo? Qual a população máxima atingida? f) Plotar num mesmo sistemas de eixos os pontos tabelados e a e o polinômio interpolador determinado no ítem b. Resposta: ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ ⎤ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ ⎡ = 65432 65432 65432 65432 65432 1313131313131 1010101010101 7777771 5555551 3333331 1111111 0000001 V Determinante Normalizado= 0.9072675023 x 10 -11 Cond (V)= 39124291.11 p(t) = .5095984263e-4*t^6-.2488977072e-2*t^5+.4187474562e-1*t^4-.2409954552*t^3- .6536635972*t^2+10.85522232*t+50. p(4) ≅ 75.91851648 população máxima ≅ p(5.312779138)= 77.05456458 41 4-Interpolação Polinomial APLICAÇÃO 2 Para determinar a resistência elétrica de um solo num sistema de aterramento, enterra-se duas hastes de cobre e aplica-se uma determinada voltagem, resultando numa corrente elétrica. Numa experiência deste tipo, foram obtidos os seguintes dados: 5,4 50 3,4 47 5,3 40 8,2 35 2 30 ))A(amperecorrente(y ))V(voltsvoltagem(x − − Estimar a corrente se a voltagem aplicada for de 43A usando o Polinômio Interpolador de Newton. Resposta: p(43)≅ 3,88 APLICAÇÃO 3 “Ao considerar que no Japão a vida média já é superior a 81 anos, a esperança de vida no Brasil de pouco mais que 71 anos ainda é relativamente baixa. E, de acordo com a projeção mais recente da mortalidade, somente por volta de 2040 o Brasil estaria alcançando o patamar de 80 anos de esperança de vida ao nascer. (Ver www.ibge.gov.br em População / Tábuas Completas de Mortalidade). A esperança de vida ao nascer de 71,3 anos coloca o Brasil na 86ª posição no ranking da ONU, considerando as estimativas para 192 países ou áreas no período 2000-2005 (World Population Prospects: The 2002 Revision; 2003). Entre 1980 e 2003 a esperança de vida ao nascer, no Brasil, elevou-se em 8,8 anos: mais 7,9 anos para os homens e mais 9,5 anos para as mulheres. Em 1980, uma pessoa que completasse 60 anos de idade teria, em média, mais 16,4 anos de vida, perfazendo 76,4 anos. Vinte e três anos mais tarde, um indivíduo na mesma situação alcançaria, em média, os 80,6 anos. Aos 60 anos de idade os diferenciais por sexo já não são tão elevados comparativamente ao momento do nascimento: em 2003, ao completar tal idade, um homem ainda viveria mais 19,1 anos, enquanto uma mulher teria pela frente mais 22,1 anos de vida”.Na tabela acima obtemos informações sobre a esperança de vida às idades exatas, especialmente, a esperança de vida ao nascer que expressa o número médio de anos que se espera viver um recém-nascido (que, ao longo da vida, estivesse exposto aos riscos de morte da tábua de mortalidade em questão 42 E.B.Hauser – Cálculo Numérico 5. Numa esfera de superfície conhecida, o coeficiente de absorção 0,7 foi mantido à temperatura de 6000 Ko . Foi calculada a energia irradiada de acordo com o tempo de irradiação, obedecendo à tabela . Pede-se para obter a possível energia irradiada quando a irradiação atingir o tempo de 25 minutos. ENERGIA IRRADIADA (Joules) TEMPO DE IRRADIAÇÃO (s) 71,72 . 310 600 94,72 . 310 800 118,4 . 310 1000 142,08 . 310 1200 165,76 . 310 1400 189,44 . 310 1600 6. Uma corda foi tensionada sob a ação de pesos distintos, quando para os respectivos pesos foram calculadas as devidas velocidades de propagação que estão indicadas abaixo. Pede-se para calcular a velocidade de propagação quando a corda está tensionada sob a ação de um peso de 7250 gf. PESO (gf) VELOCIDADE (cm/s) 6000 13728,13 6500 14288,69 7000 14828,07 7500 15348,51 8000 15851,87 Respostas 1) 126904937,1)12,0(p ≅ 2) A demanda mínima é 8632,14 Mw. e ocorre entre 3h e 4h da manhã. A demanda máxima é 101,43 Mw. e ocorre entre 14h e 15h. 3) 4128,11)850(p = 4) 4512685,13)80(p 46701783,13)80(p = = 5) 5937,177618)1500(p = 6) p(7250) = 15090,53234 45 5. Ajuste de Funções Aplicação1: Os dados acima tabelados descrevem a intensidade da luz como uma função da distância da fonte, I(d), medida num experimento. Determinar I (d) ≅ Y (d) = CBdAd 1 2 ++ . 15.0 75 18.0 70 21.0 65 24.0 60 28.0 55 34.0 50 42.0 45 52.0 40 67.0 35 85.0 30 I d Aplicação 2:Segundo a lei de Boyle, o volume de um gás é inversamente proporcional à pressão exercida, mantendo-se constante a temperatura. Para um certo gás, foram observados os seguintes valores: 04,1 5,3 15.1 0,3 28,1 5,2 43,1 0,2 60,1 5,1 91,1 0,1 27,2 5,0 85,2 0,0 Volume essãoPr Ajustar os dados tabelados a uma hipérbole do tipo: Vol(p) ≅ pB A)p(Y + = 46 E.B.Hauser – Cálculo Numérico 5.1 Introdução O ajustamento é uma técnica de aproximação. Conhecendo-se dados experimentais , , )y,x( ii , n,...,1,0i = , deseja-se obter a lei )x(fy = relacionando x com y. Devido aos erros experimentais nos n+1 pares, ))x(f,x( ii , teremos em geral 11 )x(f ε+ ; 22 )x(f ε+ ; ... ; nn )x(f ε+ , isto é , é impossível calcular exatamente a função f(x). Por isso , em vez de procurarmos a função f tal que passa por cada um dos pontos experimentais, determinaremos a função que melhor se ajusta aos pontos dados. O ajustamento traduz um comportamento médio. Para ajustar uma tabela de dados a uma função devemos: • Conhecer a natureza física do problema ; • Determinar o tipo de curva a que se ajustam os valores tabulados (graficamente e/ou cálculo das diferenças finitas ou divididas) ; • Calcular os parâmetros da curva. 5.2. Escolha da Função de Ajuste a) Função Linear (regressão linear) : xaa)x(Y 10 += , se iy∆ ≅ cte ou iy∆ ≅ cte . b) Parábola (ajuste quadrático): 2210 xaxaa)x(Y ++= , se iy 2∆ ≅ cte ou .2 cteiy ≅∆ c) Polinômio de grau p: se i p y∆ ≅ cte ou .cteyi p ≅∆ d) Função Exponencial: xab)x(Y = , se .cte x ylog i i ≅ ∆ ∆ d) FunçãoPotência: pax)x(Y = , cte xlog ylog i i ≅ ∆ ∆ Outros tipos de Funções Ajuste: • xaa y 1 xaa 1)x(Y 10 10 +=⇒ + = ; ( ) ( ) .cte x y/1y/1 i i i ≅= ∆ ∆ ∆ • ;xaa y x xaa x)x(Y 10 10 +=⇒ + = ( ) ( )( ) .ctex y/xy/x i ii ii ≅= ∆ ∆ ∆ • 47 5 - Ajuste de Funções Exemplo2 : Dada a tabela , achar a equação da reta que se ajusta usando o método dos Mínimos Quadrados. i ix iy iyix 2ix i Y ( )2iyiY − 0 0 2 0 0 2.07142857 0.005102041 1 1 3 3 1 3.14285714 0.02040816 2 2 5 10 4 4.21428571 0.61734694 3 3 5 15 9 5.2857143 0.08163266 4 4 5.5 22 16 6.3571429 0.73469755 5 5 8 40 25 7.4285714 0.32653061 ∑ 15 28.5 90 55 --------------- 1.78571440 Seja xaaxY 10)( += a função que ajusta os dados . Os parâmetros 0a e 1a constituem a solução do sistema : ∑ ∑ ∑ ∑ ∑ += ++= 2 10 10)1( ixaixaiyix ixaaniy ⇒ ⎩ ⎨ ⎧ =+ =+ 90155015 5.2811506 aa aa ⇒ ⎩ ⎨ ⎧ = = 071428571.11 071428572.20 a a Logo, xY 071428571.1071428571.2 += . 50 E.B.Hauser – Cálculo Numérico 5.3.2 - Ajuste Quadrático Seja 2210 xaxaaY ++= a função de ajuste. Pelo critério dos Mínimos Quadrados : ( ) ( )∑ ∑ −++=−= = = n 0i n 0i 2 i 2 i2i10 2 ii210 yxaxaayY)a,a,a(F assume valor mínimo se 0 a F a F a F 210 === δ δ δ δ δ δ . Assim, os parâmetros 210 aea,a são obtidos resolvendo: ( ) ⎪ ⎪ ⎪ ⎩ ⎪⎪ ⎪ ⎨ ⎧ ++= ++= +++= ∑ ∑ ∑∑ ∑ ∑ ∑ ∑ ∑ ∑∑ 4 ix2a 3 ix1a 2 ix0aiy 2 ix 3 ix2a 2 ix1aix0aiyix 2 ix2aix1a0a1niy Exemplo3 : Encontre a expressão do polinômio de 2o grau que se ajusta aos dados da tabela abaixo. i ix iy iyix 2ix iyix 2 3ix 4 ix i y∆ iy 2∆ 0 -2 -0.01 0.02 4 -0.04 -8 16 0.52 -0.21 1 -1 0.51 -0.51 1 0.51 -1 1 0.31 -0.25 2 0 0.82 0 0 0 0 0 0.06 -0.13 3 1 0.88 0.88 1 0.88 1 1 -0.07 -0.25 4 2 0.81 1.62 4 3.24 8 16 -0.32 ---------- 5 3 0.49 1.47 9 4.41 27 81 ---------- ---------- ∑ 3 3.5 3.48 19 9 27 115 ---------- ---------- ⎪ ⎩ ⎪ ⎨ ⎧ =++ =++ =++ 92115127019 48.322711903 5.32191306 aaa aaa aaa ⇒ 806285714.0201.02102142857.0)( ++−= xxxY 51 5 - Ajuste de Funções 5.3.3-Ajuste a polinômio de Grau p O ajuste a um polinômio de grau p. nppxpaxaxaaY <++++= ,... 2 210 , exige resolver o sistema linear: ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ ⎤ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ ⎡ ∑ ∑ ∑ ∑ ∑ ∑ ∑ ∑ ∑ ∑ ∑∑+ + + p2 i 1p i p i 1p i 4 i 3 i 2 ii p i 3 i 2 ii x.................xx .................................................. x...xxxx x...xxx)1n( ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ ⎤ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ ⎡ pa a a ..... ..... 1 0 = ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ ⎤ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ ⎡ ∑ ∑ ∑ iy p ix iyix iy ........ Exemplo5 – Utilizando o sistema Maple > xv:=[0,2,3,5,6,8]; > yv:=[p(0),p(2),p(3),p(5),p(6),p(8)]; > g:=fit[leastsquare[[x,y],y=a*x^3+b*x^2+c*x+d,{a,b,c,d}]]([xv,y v]); > gll :=unapply(rhs(g),x): > plots[display]({plots[pointplot]([0,p(0),2,p(2),3,p(3),5,p(5), 6,p(6),8,p(8)]),plot(gll(x), x=-.5..8.5, y=- 150..280),plot(gll(x), x=-.5..8.5, y=-150..280,thickness=2)}); 52 E.B.Hauser – Cálculo Numérico Efeito da linearização dos dadoss 5.3.4.3 -AJUSTE POR FUNÇÃO HIPERBÓLICA : xaa 1Y 10 + = . Linearizando, temos : ⎪ ⎪ ⎩ ⎪ ⎪ ⎨ ⎧ ∑=∑+∑ ∑ ∑=++ ⇒+= === = = n 0i ii n 0i 2 i1 n 0i 0i n 0i n 0i ii10 10 y/xxaax y/1xaa)1n( xaa Y 1 5.3.4.4 -OUTROS TIPOS DE FUNÇÕES DE AJUSTE 1) xaa xY 10 + = . Linearizando, temos: ⎪ ⎪ ⎩ ⎪ ⎪ ⎨ ⎧ ∑=∑+∑ ∑=∑++ ⇒+= === == n 0i i 2 i n 0i 2 i1 n 0i 0i n 0i ii n 0i i10 10 y/xxaax y/xxaa)1n( xaa Y x 2) 2 210 xaxaa 1Y ++ = . Linearizando, temos : ⎪ ⎪ ⎩ ⎪⎪ ⎨ ⎧ ∑ ∑ ∑ ∑=++ ∑ ∑ ∑=+∑+ ∑=∑+∑++ ⇒++= i 2 i 4 i2 3 i1 2 i0 ii 3 i2 2 i1i0 i 2 i2i10 2 210 y/xxaxaxa y/xxaxaxa y/1xaxaa)1n( xaxaa Y 1 3) { { 43421 2x2ax1a 2 0ay 2cxbx cxbxalnYln aeY + + ++= = ⇒ ⎪ ⎪ ⎩ ⎪⎪ ⎨ ⎧ ∑ ∑ ∑=+∑+ ∑ ∑=∑+∑+ ∑ ∑ ∑=+++ i 2 i 4 i 3 i 2 i ii 3 i 2 ii i 2 ii ylnxxcxbxaln ylnxxcxbxaln ylnxcxbaln)1n( ⇒ 0a0 aeaaln →=⇒ . 55 5 - Ajuste de Funções Exercícios: 1) Considerando: i 0 1 2 3 4 5 6 x 1 2 3 4 5 6 7 y 34 45 63 88 120 159 205 a) Mostrar que o ajuste por uma parábola é adequado. b) Ajustar os dados a uma parábola. Resposta: a) ∆2 7y ii = ∀, b) p x x x( ) , , ,= + +3 51194762 0 4047619048 30 142857142 2) Segundo o critério dos Mínimos Quadrados, qual das duas funções x558,0332,2Y1 += e x235,02 e037,2Y = melhor ajusta os dados da tabela? x -2,3 1 3,1 y 1,2 2,5 4,2 Resposta: Y2, pois 2ii1 2 ii2 )yY()yY( −<− ∑∑ 2 ii1 )yY( −∑ = 0,194 e 2ii2 )yY( −∑ = 0,0123 3) Um filme vem sendo exibido numa determinada sala de cinema por cinco semanas e a frequência semanal, (aproximada à centena mais próxima) está dada na tabela abaixo. Utilizar ajuste linear para determinar a frequência esperada na sexta semana. semana 1 2 3 4 5 frequência 5000 4500 4100 3900 3500 Resposta: Y(x) = 5280 -360x e y(6) = 3120 y = -360x + 5280 R20,9818 = 0 1000 2000 3000 4000 5000 6000 0 1 2 3 4 5 6 56 E.B.Hauser – Cálculo Numérico 4) O número de bactérias, por unidade de volume, existente em uma cultura depois de x horas é apresentado pela tabela. Ajuste os dados tabulados a uma curva exponencial da forma y =abx e avaliar y para x=7. x 0 1 2 3 4 5 6 y 32 47 65 92 132 190 275 Resposta: y x= ×32 1483 1 42694, , e y (7) = 387, 256 5) Utilizando o critério dos Mínimos Quadrados, ajustar a uma reta os dados tabulados: xi 3 5 6 8 9 11 yi 2 3 4 6 5 8 Resposta: y = -0,33328 + 0,714285x 6) A tabela abaixo fornece uma relação entre a temperatura da água e a pressão barométrica. Ajustar os dados a uma função potência. p(mm de Hg) 680 690 700 710 720 730 740 780 T( Co ) 96.92 97.32 97.71 98.11 98.49 98.88 99.26 100.73 7) A tabela abaixo fornece uma relação entre a resistência à tração do aço em função da temperatura. 1500 617 2780 485 4450 412 5260 330 5720 250 )2/( )( cmkg Cot σ Ajustar os dados a um polinômio de grau 4. y = 15,476x0,2813 R21 = 96,5 97 97,5 98 98,5 99 99,5 100 100,5 101 660 680 700 720 740 760 780 800 y = 2E-06x40,0033x - 31,8931x + 2466,48x + 47846 - R21 = 0 1000 2000 3000 4000 5000 6000 7000 0 200 400 600 80057 6 . Integração Numérica Objetivo : Calcular a ∫ b a dxxf )( , onde a função integrando )(xf ou é conhecida por sua expressão analítica ou por uma tabela de valores ( ))x(f,x ii , i = 0, 1,..., n. Aplicação: Para controlar a poluição térmica de um rio, a temperatura (oF) foi registrada, reproduzindo os dados: 1.75 17 6.78 16 1.81 15 4.86 14 5.86 13 8.84 12 1.83 11 0.77 10 3.75 9 )atemperatur(y )hora(x Encontrar a temperatura média da água entre 9h da manhã e 5h da tarde e estimar o erro cometido nesse cálculo. ]b,a[emcontínuaéf e0)b(f).a(fse],b,a[em)x(fdemédiovaloroédx)x(f ab 1fm:OBS b a >∫ − = 6.1 -Fórmulas de Newton-Cotes Consideremos ( ))x(f,x ii , i = 0, 1,..., n., hixix =−+1 , n xnxh 0−= , )( ixfiy = . A integral da função f(x) no intervalo [ ]n0 x;x é dada por : ∫ ⎟⎟ ⎠ ⎞ ⎜⎜ ⎝ ⎛ −− ++ −− + − += =∫≅∫ −n x 0x 0 n n 1n0 0 2 2 10 0 0 0 nx 0x n nx 0x dxy h!n )xx)...(xx(...y h!2 )xx)(xx(y h )xx(y dx)x(Pdx)x(f ∆∆∆ Se h xxR 0−= , então hdrdxRhxx 0Rxx 0 0 =→+= =→= e nRnxx =→= . Assim: 60 E.B.Hauser – Cálculo Numérico dR)y !n )1nR)...(1R(R...y !2 )1R(RyRy(hdR)R(Phdx)x(f 0 nnx 0x 0 2 0 n 0 0 n 0 ∆∆∆∫ ∫∫ +−− ++ − ++=≅ = ⎥ ⎥ ⎦ ⎤ ⎢ ⎢ ⎣ ⎡ +−−++−++∫ ∫ ∫∫ n 0 n 0 n 0 0 n 0 2n 000 dR)1nR)...(1R(R !n y...dR)1R(R !2 yRdRydRyh ∆∆∆ Na prática não é usual aproximar f(x) por um polinômio de grau n (elevado) devido ao erro de arredondamento que ocorre no processo. 6.2 – Regra dos Trapézios Considerando n = 1 na fórmula de Newton-Cotes temos : ∫ ∫ ∫ ∫ ⎥⎦ ⎤ ⎢⎣ ⎡ +=≅1 x 0x 1 0 1 0 1 000 RdRydRyhdR)R(Phdx)x(f ∆ [ ]1020100 ]2/R[y]R[yh ∆+= = = [ ]1000 yy2 h 2 yyh +=⎥⎦ ⎤ ⎢⎣ ⎡ + ∆ , ou para o intervalo [ ]1ii x;x + : [ ]∫ ∫+ + ++=≅ 1ix ix 1i i 1ii yy 2 hdR)R(Phdx)x(f Generalizamos para n subintervalos: [ ]n1n21on x ox y)yyy(2y 2 hdx)x(f +++++≅∫ −L Erro de Truncamento(para n subintervalos): )('' ],[ max)0(12 2 xf nxoxx xnx h TE ∈ −≤ ou iy xnx TE 2max 12 )0( ∆ − ≤ Vê-se que a fórmula dos Trapézios é exata para polinômios do 1o grau. Exemplo1: Determinar h de tal forma que a regra trapezoidal forneça o valor de dxe 1 0 2x∫ − com um erro de truncamento menor que 410− . 0245.0h10)2( 12 h)x(''fmáx)01( 12 hE 4 22 T <⇒<=−≤ − [ ]1;0x∈ h/)xx(n 0n −= , n/)xx(h 0n −= , 0245.0n/)xx( 0n <− , 41n8.40n =→> , 02439.041 1h == 61 6 -Integração Numérica 6.3 – Fórmula de Simpson Fazendo n = 2 na fórmula de Newton-Cotes, temos, [ ]∫ ∫ ∫∫ ++= ⎥ ⎥ ⎦ ⎤ ⎢ ⎢ ⎣ ⎡ −++≅2 x 0x 210 2 0 2 0 0 22 000 yy4y 3 hdR)1R(R !2 yRdRydRyhdx)x(f ∆∆ Generalizamos para n subintervalos, n par: [ ]∫ +++++++++++≅ −−nx ox n2n6421n531o y)yyyy(2)yyyy(4y3 hdx)x(f LL ERRO DE TRUNCAMENTO PARA A FÓRMULA DE SIMPSON )x(fmax)xx( 180 hE '''' ]nx,ox[x 0n 4 S ∈ −≤ ou i 40n S ymax180 )xx(E ∆−≤ Exercícios 1.Aplicação:.Para controlar a poluição térmica de um rio, a temperatura (oF) foi registrada, reproduzindo os dados: 1.75 17 6.78 16 1.81 15 4.86 14 5.86 13 8.84 12 1.83 11 0.77 10 3.75 9 )atemperatur(y )hora(x Encontrar a temperatura média da água entre 9h da manhã e 5h da tarde e estimar o erro cometido nesse cálculo. 2. Estimar a área da região hachurada pela regra dos Trapézios e pela de Simpson, usando oito subintervalos. 62 7 - Resolução Numérica de Equações Diferenciais Ordinárias Objetivo: Resolver numericamente(e generalizar para problemas de ordem mais elevada) o problema de valor inicial de primeira ordem (PVI) ⎪⎩ ⎪ ⎨ ⎧ = = oo y)x(y )y,x(f'y , e o problema de valor de contorno de segunda ordem, linear: (PVC) ⎪⎩ ⎪ ⎨ ⎧ == =++ .y)x(y,y)x(y )x(fy)x(Q'y)x(P''y nnoo Os métodos numéricos são processos que fornecem valores aproximados da solução em pontos particulares, usando operações algébricas. Os métodos que estudaremos determinarão estimativas da solução nos pontos K,x,x,x 21o , onde ihxx oi += , i= 0,1,...,n . A escolha do valor de h é arbitrária e, em geral, quanto menor h , melhor a estimativa da solução obtida. Sejam iyixy ≅)( , ihxx oi += , i= 0,1,...,n e n/)xx(h on −= . Método de Euler(PVI): )y,x(hfyy iii1i +=+ , i= 0,1,...,n-1. Método de Runge-Kutta de 2 a ordem(PVI) : ( )21i1i kk2 hyy ++=+ , )y,x(fk ii1 = e )hky,hx(fk 1ii2 ++= , i= 0,1,...,n-1. Diferenças Finitas(p/diferenças centrais)(PVC): Para i= 1,...,n-1 , é gerado um sistema de n-1 equações lineares: ( ) )x(fhy)x(P 2 h1y)x(Qh2y)x(P 2 h1 i 2 1iiii 2 1ii =⎟⎠ ⎞ ⎜ ⎝ ⎛ −++−+⎟ ⎠ ⎞ ⎜ ⎝ ⎛ + −+ Sistemas de Equações Diferenciais de primeira ordem com condições iniciais ⎪ ⎪ ⎩ ⎪⎪ ⎨ ⎧ == = = oooo z)x(z,y)x(y )z,y,x(g'z )z,y,x(f'y Para obtermos sua solução é possível aplicar os métodos de Euler e Runge-Kutta de segunda ordem. Por exemplo, por Euler as estimativas são obtidas aplicando. yi+1 = yi + h f (xi, yi, zi ) e zi+1 = zi + h g (xi, yi, zi ). Mudança de Variável para problemas de valor inicial de segunda ordem: Para ⎪⎩ ⎪ ⎨ ⎧ == = oooo z)x('y,y)x(y )z,y,x(f''y , fazemos a mudança de varável u'y = , com oo y)x(y = . Então devemos resolver o sistema de duas equações diferenciais ordinárias, acopladas pelas condições iniciais: ⎪ ⎪ ⎩ ⎪⎪ ⎨ ⎧ == = = oooo z)x(u,y)x(y )z,y,x(f'u u'y . 65 E.B.Hauser – Cálculo Numérico Derivada Diferença Finita h, k=tamanho do passo na direção x e na direção y (ou t) )ix´(y atrasada h 1iyiy,centrada h2 1iy1iy,avançada h iy1iy −−−−+−+ )ix´´(y 2h 1iyiy21iy −+−+ centrada )ix('''y 3h2 2iy1iy21iy22iy −−−++−+ centrada )ix( IVy 4h 2iy1iy4iy61iy42iy −+−−++−+ centrada ( )jt,ixt u ∂ ∂ k2 1j,i u 1j,i u − − + centrada ( )jt,ix2x u2 ∂ ∂ 2h j,1i u j,i u2 j,1i u − +− + centrada ( )jy,ix2y u2 ∂ ∂ 2k 1j,i u j,i u2 1j,i u − +− + centrada Exemplos: 66 7 -Resolução Numérica de Equações Diferenciais Ordinárias Método de Euler Problema: y' = y – x; y(0) = 2 yn xn h = 0.1 h = 0.05 h = 0.01 h = 0.005 Solução exata Y(x) = ex + x + 1 0.0 2.0000 2.0000 2.0000 2.0000 2.0000 0.1 2.2000 2.2025 2.2046 2.2049 2.2052 0.2 2.4100 2.4155 2.4202 2.4208 2.4214 0.3 2.6310 2.6401 2.6478 2.6489 2.6499 0.4 2.8641 2.8775 2.8889 2.8903 2.8918 0.5 3.1105 3.1289 3.1446 3.1467 3.1487 0.6 3.3716 3.3959 3.4167 3.4194 3.4221 0.7 3.6487 3.6799 3.7068 3.7102 3.7138 0.8 3.9436 3.9829 4.0167 4.0211 4.0255 0.9 4.2579 4.3066 4.3486 4.3541 4.3596 1.0 4.5937 4.6533 4.7048 4.7115 4.7183 Método de Runge-Kutta de segunda ordem Problema: y' = y – x; y(0) = 2 yn xn h = 0.1 h = 0.05 h = 0.01 Solução exata Y(x) = ex + x + 1 0.0 2.0000000 2.0000000 2.0000000 2.0000000 0.1 2.2050000 2.2051266 2.2051691 2.2051709 0.2 2.4210250 2.4213047 2.4213987 2.4214028 0.3 2.6492326 2.6496963 2.6498521 2.6498588 0.4 2.8909021 2.8915852 2.8918148 2.8918247 0.5 3.1474468 3.1483904 3.1487076 3.1487213 0.6 3.4204287 3.4216801 3.4221007 3.4221188 0.7 3.7115737 3.7131870 3.7137294 3.7137527 0.8 4.0227889 4.0248265 4.0255115 4.0255409 0.9 4.3561818 4.3587148 4.3595665 4.3596031 1.0 4.7140809 4.7171911 4.7182369 4.7182818 67 E.B.Hauser – Cálculo Numérico Exercícios: 1. Utilizando o método de Euler, determinar y(X i ), se : a) ⎪⎩ ⎪ ⎨ ⎧ = −−= 2)1(y 1xyỳ , h=0,25 , Xi =1,75 b) ⎪⎩ ⎪ ⎨ ⎧ = >+= 0)1(y 0y,yxỳ , h=0,1 , Xi =1,4 c) ⎪⎩ ⎪ ⎨ ⎧ = += 2,0)1(y y2xỳ 2 , h=0,1 , Xi =1,4 d) ⎪⎩ ⎪ ⎨ ⎧ = −−= 2)0(y 100e101y100ỳ x , h=0,1 , Xi =0,3 e) ⎪⎩ ⎪ ⎨ ⎧ −== =−− 5,0)1`(y,5.0)1(y 0yx8ỳ``xy 33 , h=0,1 , Xi =1,3 f) ⎪⎩ ⎪ ⎨ ⎧ == =+−− 0)0`(y,5,0)0(y 0yỳ)y1(``y 2 , h=0,1 , Xi =0,3 2.Utilizando o Método de Heun (Runge-Kutta de 2ª ordem), determinar y( Xi ), se: a) ⎪⎩ ⎪ ⎨ ⎧ = += 1)0(y y2x3ý , h=0,2, Xi =0,4 b)⎩ ⎨ ⎧ = −= 1)0(y xyý , h=0,1, Xi =0,2 c) ⎪⎩ ⎪ ⎨ ⎧ = = 0)0(y 3x4ý , h=0,1 , Xi =0,3 3.Utilizando o Método das diferenças finitas, e o valor indicado de n, resolver o PVC. a) 4n, 1)2(y,4)0(y 0y9''y = ⎩ ⎨ ⎧ == =+ b) 5n, 0)1(y,4)0(y x5y'y2''y = ⎩ ⎨ ⎧ == =++ c) 8n, 0)2(y,5)1(y 0y3'xy3''yx2 = ⎪⎩ ⎪ ⎨ ⎧ == =++ d) 10n, 2)1(y,0)0(y xxy'y)x1(''y = ⎩ ⎨ ⎧ == =+−+ 70 7 -Resolução Numérica de Equações Diferenciais Ordinárias 4. PVI -Considerar um sistema massa-mola-amortecedor descrito pela equação diferencial ordinária de segunda ordem: m y” + cy’ + ky = L sin x. Utilizando Euler com h=0.25, estimar o deslocamento para o tempo x=0.5, para massa m=1, amortecedor c=0.5, rigidez k=2, amplitude da força L=0.5 , com deslocamento inicial y(0)=1, velocidade inicial y’(0) =0. 5. PVC - Considerar o problema de deflexão de uma viga de seção transversal retangular sujeita a uma carga uniforme, tendo seus extremos apoiados de modo a não sofrer deflexão alguma. O problema de valor de contorno que rege essa situação física é Lx0),1x(x EI2 qw EI S dx wd 2 2 <<−+= , Como não ocorre deflexão nas extremidades da viga, as condições de contorno são w(0) =0, w(L) =0. Considerando: • Comprimento L=120 pol; • Intensidade de carga uniforme q=100 lb/pé; • Módulo de elasticidade E=3.0 x 107 lb/pol2; • Esforço nas extremidades S=1000 lb; • Momento central de Inércia I=625 pol4; aproximar a deflexão w(x) da viga a cada 20pol, utilizando diferenças finitas. Respostas : 1.a) y(1,75)=1,2018 2.a) y(0,4)=1,1189 b) y(1,4)=0,4558 b) y(0,2)=0,9801 c) y(1,4)=0,7778 c) y(0,3)=0,009 d) y(0,3)=-25,2206 e) y(1,3)= 0, 3647 f) y(0,3)= 0,3991 3.a) 3226,6y 5807,2y 6774,5y 3 2 1 = −= −= b) 2167,0y 3308,0y 3356,0y 2259,0y 4 3 2 1 −= −= −= −= c) 2913,0y 6430,0y 0681,1y 5826,1y 2064,2y 9640,2y 8842,3y 7 6 5 4 3 2 1 = = = = = = = d) 8474,1y 6855,1y 5149,1y 3353,1y 1465,1y 9471,0y 7357,0y 5097,0y 2660,0y 9 8 7 6 5 4 3 2 1 = = = = = = = = = 4) y(0.5)=0.875 71 8 – Resolução Numérica de Equações Diferenciais Parciais 8.1 -Introdução Equação diferencial parcial (EDP) é a uma equação que envolve duas ou mais variáveis independentes ),t,z,y,x( K e derivadas parciais de uma função incógnita(variável dependente que queremos determinar) ),t,z,y,x(uu K≡ . Um corpo é isotrópico se a condutividade térmica em cada um de seus pontos é independente da direção do fluxo de calor através do ponto. Em um corpo isotrópico, a temperatura, ),,,,( tzyxuu ≡ é obtida resolvendo-se a equação diferencial parcial (EDP) t ucp z uk zy uk yx uk x ∂ ∂ =⎟ ⎠ ⎞ ⎜ ⎝ ⎛ ∂ ∂ ∂ ∂ +⎟⎟ ⎠ ⎞ ⎜⎜ ⎝ ⎛ ∂ ∂ ∂ ∂ +⎟ ⎠ ⎞ ⎜ ⎝ ⎛ ∂ ∂ ∂ ∂ onde k, c e p são funções de (x,y,z), e representam respectivamente, a condutividade térmica, o calor específico e a densidade do corpo no ponto (x,y,z). Quando k, c e p são constantes, essa equação é denominada equação simples tridimensional do calor, e é expressa como . t u k cp z u y u x u 2 2 2 2 2 2 ∂ ∂ = ∂ ∂ + ∂ ∂ + ∂ ∂ Se o domínio do problema é relativamente simples, a solução dessa equação é obtida utilizando a série de Fourier. Na maioria das situações onde k, c e p não são constantes ou quando o domínio é irregular, a solução da equação diferencial parcial deve ser obtida por meio de métodos de aproximação. Para introduzir métodos numéricos de resolução de EDP, utilizaremos as equações de Poisson, do Calor e da Onda, as quais representam protótipos das EDP´s elípticas, parabólicas e hiperbólicas. Será adotado um procedimento geral, seguindo os passos: 1) Construir uma malha a partir do domínio do problema; 2) Para os pontos interiores da malha, escolher a discretização das derivadas parciais; 3) Construir o sistema de equações lineares usando a discretização dos pontos interiores, ( )ji y,xf e as condicções de contorno. 4) Resolver o sistema de equações lineares(escolher o método masi eficiente), cuja solução forne as aproximações da solução nos pontos interiores da malha. 8.1.1-Equação Do Potencial ou de Poisson(EDP Elíptica) Consideremos a equação de Poisson: ).y,x(f)y,x( y u)y,x( x u 2 2 2 2 = ∂ ∂ + ∂ ∂ 72 E.B.Hauser – Cálculo Numérico )()0,( xfxu = e ),x(g)0,x( t u = ∂ ∂ Se os pontos extremos forem fixos, teremos: 0)t,l(ue0)t,0(u == . Os outros problemas físicos envolvendo a equação diferencial parcial hiperbólica ocorrem no estudo de vigas vibrantes com uma ou ambas as extremidades clamped e na transmissão de eletricidade em uma linha de transmissão longa onde exista alguma perda de corrente para o solo. 8.2. -Método das Diferenças Finitas para Equação Diferencial Parcial Elíptica Consideremos o problema de valor de contorno envolvendo a equação de Poisson (que é uma equação diferencial parcial elíptica) ⎪ ⎩ ⎪ ⎨ ⎧ ∈= = ∂ ∂ + ∂ ∂ ≡∇ ,S)y,x(para)y,x(g)y,x(u )y,x(f)y,x( y u)y,x( x u)y,x(u 2 2 2 2 2 onde S é a fronteira(contorno) do retângulo }dyc,bxa/)y,x{(R <<<<= . Se f e g são contínuas em seus domínios, então existe uma única solução para esse problema de valor de contorno.. Utilizaremos uma adaptação do método de Diferenças Finitas . Em R contruímos uma grade(figura4) , traçando linhas verticais e horizontais ( ixx = e , jyy = grid lines) pelos pontos ( ),, ji yx e suas intersecções são os pontos de rede (mesh points), onde ,ihaxi += h=(b – a)/n , para todo i = 0,1,...,n e ,jkcy j += k=(d – c)/m , para todo j = 0,1,...,m Figura 4 75 8 - Resolução Numérica de Equações Diferenciais Parciais Para cada ponto de rede no interior da quadricula ( ),, ji yx com i = 1,2,...,n-1 e com j = 1,2,...,m-1, utilizamos a série Taylor na variável x ao redor de ix para gerar a fórmula das diferenças centrais ( ) ( ) ( ) ( ) ( ),y, x u 12 h h y,xy,xu2y,xu y,x x u ji4 42 2 j1ijij1i ji2 2 ξ ∂ ∂ − +− = ∂ ∂ −+ (2.1) onde ( )., 11 +−∈ ii xxξ Também utilizamos a série de Taylor na variável y ao redor de jy para gerar a fórmula das diferenças centrais ( ) ( ) ( ) ( ) ( ),,x y u 12 k k y,xy,xu2y,xu y,x y u ji4 42 2 j1ijij1i ji2 2 η ∂ ∂ − +− = ∂ ∂ −+ (2.3) onde ( )., 11 +−∈ jj yyη A utilização dessas formulas na Equação (2.1) nos permite expressar a equação de Poisson nos pontos ( ),, ji yx como ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ),,x y u 12 ky, x u 12 hy,xf k y,xy,xu2y,xu h y,xy,xu2y,xu ji4 42 ji4 42 ji 2 j1ijij1i 2 j1ijij1i ηξ ∂ ∂ + ∂ ∂ += +− + +− −+−+ para todo i = 1,2,...,n-1 e j = 1,2,...,m-1, e as condições de limite como ),(),( 00 jj yxgyxu = e ),,(),( jnjn yxgyxu = para todo j = 0,1,...,m; ),(),( 00 yxgyxu ii = e ),,(),( mimi yxgyxu = para todo i = 1,2,...,n-1; Na forma da equação de diferenças, isso resulta no método das Diferenças Finitas para a equação de Poisson, com um erro local de truncamnto da ordem de ( ):22 khO + ( ) ( ) ( ) 1.-m1,2,..., j e 1-n1,2,..., i ,y,xfhww k hwww1 k h2 ji 2 1j,i1j,i 2 j,1ij,1iij 2 == −=+⎟ ⎠ ⎞ ⎜ ⎝ ⎛−+− ⎥ ⎥ ⎦ ⎤ ⎢ ⎢ ⎣ ⎡ +⎟ ⎠ ⎞ ⎜ ⎝ ⎛ −+−+ (2.4) ),( 00 jj yxgw = e ),( jnnj yxgw = , para todo j = 0,1,...,m; (2.5) ),( 00 yxgw ii = e ),,( miim yxgw = para todo i = 1,2,...,n-1; onde ijw é uma aproximação de ),( ji yxu . 76 E.B.Hauser – Cálculo Numérico A equação em (2.4) envolve aproximações a ),( yxu nos pontos ( )ji yx ,1− , ( )ji yx , , ( )ji yx ,1+ , ( )1, −ji yx e ( )1, +ji yx . Reproduzindo a parte da malha na qual esses pontos estão situados (veja figura.5), observamos que cada equação contém aproximações em uma região em forma de estrela ao redor de ( )ji yx , . Figura 5 Se utilizarmos a informação das condições de limite (2.5) sempre que for conveniente no sistema dado por (12.4), poderemos dizer que, em todos os pontos ( )ji yx , adjacentes ao ponto de rede do limite, teremos um sistema linear ( n – 1)(m – 1) x ( n – 1)(m – 1), cujas incógnitas são as aproximações ijw a ),( ji yxu no interior dos pontos de rede. O sistema linear que contém essas incógnitas será expresso mais eficientemente em cálculos matriciais se for introduzida uma remarcação dos pontos interiores da malha. Um sistema de marcação desses pontos consiste em utilizar ( )jil y,xP = e ijl ww = onde ),1)(1( −−−+= njmil para todo i = 1,2,...,n-1 e j = 1,2,...,m-1. Isso marca consecutivamente os ponto de rede da esquerda para direita e de cima para baixo. Por exemplo, com n = 4 e m = 5, com a remarcação se obtém uma quadrícula cujos pontos são mostrados na Figura 12.6.Ao marcar os pontos desse modo, se garante que o sistema necessário para determinar ijw seja uma matiz de banda com uma largura de banda de, no máximo, 2n – 1. Figura 6 77 8 - Resolução Numérica de Equações Diferenciais Parciais > A := matrix( [[4,-1,0,-1,0,0,0,0,0],[-1,4,-1,0,-1,0,0,0,0],[0,- 1,4,0,0,-1,0,0,0],[-1,0,0,4,-1,0,-1,0,0],[0,-1,0,-1,4,-1,0,- 1,0],[0,0,-1,0,-1,4,0,0,-1],[0,0,0,-1,0,0,4,-1,0],[0,0,0,0,-1,0,- 1,4,-1],[0,0,0,0,0,-1,0,-1,4]]): > b := vector( [25,50,150,0,0,50,0,0,25]): > linalg[linsolve](A, b); > W:=evalf(%); Tabelamos os resultados: i 1 2 3 4 5 6 7 8 9 iw 18.75 37.50 56.25 12.5 25.00 37.50 6.25 12.50 18.75 Assim, para cada i, )( ii Puw = , isto é, iw é uma estimativa da solução em )y,x(P iii = . Por exemplo, 5.37)25.0,375.0(u)P(uw 66 ≅== Exemplo2: Consideremos o modelo para determinar a temperatura de estado estacionário para uma placa quadrada com condições de contorno dadas: ⎪ ⎪ ⎪ ⎩ ⎪ ⎪ ⎪ ⎨ ⎧ ⎩ ⎨ ⎧ <≤− << == −== <<<<= ∂ ∂ + ∂ ∂ 2x1,x2 1x0,x )2,x(u,0)0,x(u )y2(y)y,2(u,0)y,0(u 2y0,2x0,0)y,x( y u)y,x( x u 2 2 2 2 80 E.B.Hauser – Cálculo Numérico Utilizando diferenças finitas centrais com h=k=2/3, para i, j = 0,1,2,3, obtemos a equação de diferenças finitas: 0u4uuuu j,i1j,ij,1i1j,ij,1i =−+++ −−=+ Para determinar o valor de 22u)3/4,3/4(ue21u)3/2,3/4(u,12u)3/4,3/2(u,11u)3/2,3/2(u ≅≅≅≅ Utilizando o sistema Maple para resolver o sistema linear tridiagonal: ⎪ ⎪ ⎩ ⎪ ⎪ ⎨ ⎧ −=−+ −=+− −=+− =++− 9/1422u412u21u 3/222u12u411u 9/822u21u411u 012u21u11u4 > solve({-4*u11+u21+u12=0, u11-4*u21+u22=-8/9, u11-4*u12+u22=-2/3, u21+u12-4*u22=-14/9},{u11,u12,u21,u22}); { }, , , = u12 13 36 = u11 7 36 = u22 7 12 = u21 5 12 > evalf(%); { }, , , = u12 .3611111111 = u11 .1944444444 = u22 .5833333333 = u21 .4166666667 Exercícios: 1) Com h=k=10, resolver numericamente equação de Laplace ,60y0,80x0,0)t,x(2x u2)t,x(2y u2 <<<<= ∂ ∂ + ∂ ∂ sujeita às condições 100)y,0(ue,0)y,80(u)60,x(u)0,x(u ==== . j↓ 0 1 2 3 4 5 6 7 8 i← 6 100 0 0 0 0 0 0 0 0 60 5 100 46,993 24,835 13,978 8,068 4,633 2,523 1,110 0 50 4 100 63,138 38,368 23,010 13,660 7,942 4,348 1,917 0 40 3 100 67,192 42,490 26,032 15,620 9,128 5,009 2,211 0 30 2 100 63,138 38,368 23,010 13,660 7,942 4,348 1,917 0 20 1 100 46,993 24,835 13,978 8,068 4,633 2,523 1,110 0 10 0 100 0 0 0 0 0 0 0 0 0 x → 0 10 20 30 40 50 60 70 80 y↓ 2) Resolver numericamente equação de calor 4 ij,iuij,iuj,1iuj,1iu j,iu ++−+−++= 81 8 - Resolução Numérica de Equações Diferenciais Parciais ,t0,1x0,0)t,x( x u)t,x( t u 2 2 ≤<<= ∂ ∂ − ∂ ∂ com as condições de contorno ,t0,0)t,1(u)t,0(u <== e as condições iniciais .1x0),x(sen)0,x(u ≤≤= π Ajuda: Diferenças finitas para o problema com h=0,2, k=0.01 ⎟ ⎠ ⎞ ⎜ ⎝ ⎛ − + + +=+ j,1iuj,1iu25,0j,iu5,01j,iu j ↓ 0 1 2 3 4 5 ← i 6 0 0,3219 0,5208 0,5208 0,3219 0 0,06 5 0 0,3559 0,5758 0,5758 0,3559 0 0,05 4 0 0,3934 0,6366 0,6366 0,3934 0 0,04 3 0 0,4350 0,7038 0,7038 0,4350 0 0,03 2 0 0,4809 0,7781 0,7781 0,4809 0 0,02 1 0 0,5317 0,8602 0,8602 0,5317 0 0,01 0 0 0,5878 0,9511 0,9511 0,5878 0 0,00 x → 0 0,2 0,4 0,6 0,8 1 ↑ tempo Diferenças finitas para o problema com h=0,1, k=0.005 ⎟ ⎠ ⎞ ⎜ ⎝ ⎛ − + + +=+ j,1iuj,1iu1,0j,iu9,01j,iu j ↓ 0 1 2 3 4 5 6 7 8 9 10 ← i 6 0 0 0,5189 0,9869 1,3584 1,5969 1,6791 1,5969 1,3584 0,9869 0,5189 0 0,030 5 1 0 0,4759 0,9053 1,2460 1,4647 1,5401 1,4647 1,2460 0,9053 0,4759 0 0,025 4 2 0 0,4365 0,8304 1,1429 1,3435 1,4127 1,3435 1,1429 0,8304 0,4365 0 0,020 3 3 0 0,4004 0,7616 1,0483 1,2324 1,2958 1,2324 1,0483 0,7616 0,4004 0 0,015 2 4 0 0,3673 0,6986 0,9616 1,1304 1,1886 1,1304 0,9616 0,6986 0,3673 0 0,010 1 5 0 0,3369 0,6408 0,8820 1,0369 1,0902 1,0369 0,8820 0,6408 0,3369 0 0,005 0 6 0 0,3090 0,5878 0,8090 0,9511 1,0000 0,9511 0,8090 0,5878 0,3090 0 0,000 x → x 0 0,1 0,2 0,3 0,4 0,5 0,6 0,7 0,8 0,9 1 ↑ tempo 8.4 – Bibliografia para EDP`s 82 Bibliografia ATKINSON, K. E. An Introduction to Numerical Analysis. 2.ed. New York : John Wiley & Sons, 1989. AYYUB, Bilal M.; McCUEN, Richard H. Numerical Methods For Engineers. New Jersey: Prentice Hall, 1996. BARROSO, Leonidas Conceição et al. Cálculo Numérico com Aplicações. 2.ed. São Paulo : Harbra, 1987. BURDEN, Richard L., FAIRES, J. Douglas. Análise Numérica. São Paulo : Thomson, 2003. CHAPRA, Steven C., CANALE, Raymond P. Numerical Methods for Engineers with Programming and Software Applications. 4.ed. Boston : McGraw-Hill, 2002. CLÁUDIO, Dalcidio Moraes, MARINS, Jussara Maria. Cálculo Numérico Computacional. 2.ed. São Paulo : Atlas, 1994. CUNHA, M.Cristina, Métodos Numéricos, Campinas, SP, 2003, Ed. Unicamp. GANDER, Walter.; HREBICEK, Jiri. Solving Problems in Scientific Computing Using Maple and Matlab: Berlin, Springer-Verlag, 1995. HUMES, Ana Flora P. de Castro et al. Noções de Cálculo Numérico. São Paulo: McGraw- Hill, 1984. KREYSZIG, Erwin. Advanced Engineering Mathematics. New York, NY : John Wiley & Sons, 1993. O’NEIL,PeterV. Advanced Engineering Mathematics. 4.ed. Pacific Grove, CA : Brooks/Cole, 1995. RICE, Richard G.; DO, Duong D. Applied Mathematics and Modelling for Chemical Engineers. New York : John Wiley & Sons, 1995. RUGGIERO, Márcia A. Gomes, LOPES, Vera Lúcia da Rocha. Cálculo numérico: aspectos teóricos e computacionais. São Paulo : McGraw-Hill, 1997. SCHELEIDER, Maria Amélia N, Cunha, Maria Cristina dC, Métodos Numéricos para Equações Diferenciais Parciais, Notas em Matemática Aplicada, Volume 4, São Carlos, SP, 2003, SBMAC. (Disponível em http://www.sbmac.org.br/boletim/pdf_2003/livro_04_2003.pdf) STROUD, K.A, BOOTH, Dexter J., Advanced Engineering Mathematics, New York, Palgrave Macmillan, 2003. 85 Cálculo Numérico - Formulário 1a) # F = 1)1mM(1t)1(2 +⎥⎦ ⎤ ⎢⎣ ⎡ +−−− ββ , ),,,( MmtF β ; 1b) ⎥ ⎥ ⎦ ⎤ ⎢ ⎢ ⎣ ⎡ ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ + −+++−=+ 1kx kx1kxlog3,0)1kx,kx(DIGSE µ 2) Seja 0ax1a 2x2a 3x3a... 1nx1na nxna)x(p +++++ − −+= a) Horner: 0ax)1ax)2a...x)2nax)1naxna 1n (((...()x(p +++−+−+ − = 321 b) Huat: Se p(0) ≠ 0 e para algum k = 1, 2, ...,n-1, 2ka ≤ 1k1k aa +− , então p tem raízes complexas. 3) Método de Newton-Raphson: 0)kx('fmáx,...,2,1,0k,)kx('f )kx(f kx1kx ≠=−=+ 4) Sistema de n equações lineares: AX=B, A = ⎥⎦ ⎤ ⎢⎣ ⎡ ija , X = ⎥⎦ ⎤ ⎢⎣ ⎡ ix e B = ⎥⎦ ⎤ ⎢⎣ ⎡ ib com i,j = 1, 2, ..., n. 4a) 2kna 2 2ka 2 1kakonden21 AdetANORM L L ++== α ααα , para k = 1, 2, ..., n. 4b) A matriz A é Diagonal Dominante se .n,..,2,1j,i n ij 1j ijaiia =∀ ≠ = > ∑ 4c)Gauss -Jacobi e Gauss -Seidel: ∀ i= 1, 2, ..., n e k = 1, 2, ..., máx, ⎟ ⎟ ⎟ ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎜ ⎜ ⎜ ⎝ ⎛ += − + − = −= + ⎟ ⎟ ⎟ ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎜ ⎜ ⎜ ⎝ ⎛ ≠ = −= + ∑∑∑ kjx n 1ij ija1kj xija 1i 1j ib iia 1 1ki x kj xija n ij 1j ib iia 1 1ki x 5) Diferenças Finitas Ascendentes: ∀ i= 0, 1, 2, ..., n ,seja ii 0 yy =∆ . Para k = 1, 2, ..., n a diferença finita de ordem k é iy 1k 1iy 1k iy k −−+ −= ∆∆∆ , com i= 0, 1, 2, ..., n-k . 6) Diferenças Divididas: ∀ i= 0, 1, 2, ..., n ,seja ii 0 yy =∆ . Para k = 1, 2, ..., n, a diferença dividida de ordem k é ixkix iy 1k 1iy 1k iy k −+ −−+ − = ∆∆ ∆ , com i= 0, 1, 2, ..., n-k . 7) Polinômio Interpolador de Newton para Diferenças Finitas Ascendentes: oy n nh!n )1nxx()1xx)(oxx( oy 2 2h!2 )1xx)(oxx( 0yh )oxx( oy)x(p ∆∆∆ −−−−++ −− + − += L L Se h )oxx(z − = , !n oy n ))1n(z()1z(z !3 oy 3 )2z)(1z(z !2 oy 2 )1z(z0yzoy)z(p ∆∆∆ ∆ −−−++−−+−++= LL . 8) Polinômio Interpolador de Newton para Diferenças Divididas oy n)1nxx)...(1xx)(oxx(...oy 2)1xx)(oxx(0y)oxx(oy)x(p ∆∆∆ −−−−++−−+−+= 86 Cálculo Numérico - Formulário 9) Ajustamento a um Polinômio de grau p: Se pxpa 2x2ax1aoa)x(Y ++++= L é a função que ajusta os pontos npen,...1,0i,)y,x( ii <= , então, para ∑ ∑ = = n 0i , os parâmetros pa,...,1a,0a constituem a solução de: ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎝ ⎛ = ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎝ ⎛ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎝ ⎛ +++ + + ∑ ∑ ∑ ∑∑∑∑∑ ∑∑∑∑∑ ∑∑∑∑ iy p ix iyix iy pa 1a oa p2 ix 3p ix 2p ix 1p ix p ix 1p ix 4 ix 3 ix 2 ixix p ix 3 ix 2 ixix1n MM L MLMMMM L L 9a)Ajuste Linear e Quadrático x1aoa)x(Y += , e 2x2ax1aoa)x(Y ++= ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ ⎟⎟ ⎟ ⎠ ⎞ ⎜⎜ ⎜ ⎝ ⎛ + ∑∑ ∑ 1a 0a 2 ixix ix1n = ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ ∑ ∑ iyix iy , e ⎟ ⎟ ⎟ ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎜ ⎜ ⎜ ⎝ ⎛ + ∑∑∑ ∑∑∑ ∑∑ 4 ix 3 ix 2 ix 3 ix 2 ixix 2 ixix1n ⎟⎟ ⎟ ⎟ ⎟ ⎠ ⎞ ⎜⎜ ⎜ ⎜ ⎜ ⎝ ⎛ 2a 1a oa = ⎟ ⎟ ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎜ ⎜ ⎝ ⎛ ∑ ∑ ∑ iy 2 ix iyix iy 10) Integração: Para noxnxh /)( −= , ihoxix += e )( ixfiy = , i = 0,1,...n, 10a) Trapézios: [ ]n1n21o x x y)yyy(2y 2 hn o dx)x(f +++++≅ −∫ L )x(''f no max)xx( 12 hE ]x,x[x 0n 2 T ∈ −≤ ou i 20n T ymax12 )xx(E ∆−≈ . 10b) Simpson ( n par) : ∫ ⎥⎦ ⎤ ⎢⎣ ⎡ +−+++++−+++++≅ nx ox ny)2ny6y4y2y(2)1ny5y3y1y(4oy3 hdx)x(f LL )x(''''f ]nx,ox[x max)0xnx(180 4h SE ∈ −≤ ou i 40n S ymax180 )xx(E ∆−≈ . 11) PVI: Sejam ⎪⎩ ⎪ ⎨ ⎧ = = oy)ox(y )y,x(f'y iyixy ≅)( com ihoxix += , i= 0,1,...,n e n/)oxnx(h −= . 11 a)Euler: )iy,ix(hfiy1iy +=+ , i= 0,1,...,n-1. 11b) Runge-Kutta de 2 a Ordem : ( )2k1k2 h iy1iy ++=+ , )iy,ix(f1k = e )hky,hx(fk 1ii2 ++= , i= 0,1,...,n-1. 12) PVC: Sejam ⎪⎩ ⎪ ⎨ ⎧ == =++ .y)x(y,y)x(y )x(fy)x(Q'y)x(P''y nnoo , iyixy ≅)( com ihoxix += , i= 0,1,...,n e n/)oxnx(h −= . Diferenças Finitas: ( ) )x(fhy)x(P 2 h1y)x(Qh2y)x(P 2 h1 i 2 1iiii 2 1ii =⎟⎠ ⎞ ⎜ ⎝ ⎛ −++−+⎟ ⎠ ⎞ ⎜ ⎝ ⎛ + −+ , i= 1,...,n-1 87 Aula de Laboratório 1 utilizando Maple Operações, Funções Básicas, Constantes Notação Exemplos Adição + > 3+5; Subtração - > 758-195; Multiplicação * > 7.5*8; Divisão / > 39/13; Potenciação ^ > 2^12; Valor Absoluto de x abs ( x ) > abs(-7); Raiz Quadrada de x sqrt (x) >sqrt(7835); > sqrt(7835.); Raiz n-ésima de x x^(1/n) > 27^(1/3); Fatorial de n n ! > 28!; π Pi > Pi; Infinito infinity > infinity; Unidade Imaginária sqrt(-1) ou I > sqrt (-1); Número de Euler: e exp(1) ou E > exp(1); Função Exponencial exp(x) b^x >exp(3.5); >7^x; Logaritmo Natural Logaritmo de base b ln(x) log[b](x) = ln(x)/ln(b). > ln (7); > log[3](10); Funções Trigonométricas , Hiperbólicas e suas inversas (argumento x em Radianos) sin(x) cos(x) tan(x) sec(x) csc(x) cot(x) sinh(x) cosh(x) tanh(x) sech(x) csch(x) coth(x) arcsin(x) arccos(x) arctan(x) arcsec(x) arccsc(x) arccot(x) arcsinh(x) arccosh(x) arctanh(x) arcsech(x) arccsch(x) arccoth(x) > sin(Pi/2); > cos(75.3); >arccosh(0); Alguns Comandos 1) > evalf(expr); avalia expr utilizando aritmética de ponto flutuante com precisão determinada pela variável global Digits. > evalf(27^(1/3)); > evalf( Pi); > evalf(arccosh(0)); 2) > evalf(expr, n); calcula expr com n dígitos de precisão. > evalf( Pi, 21); > evalf(ln (7),8); 3) > Digits := n; ajusta para n o número de dígitos utilizados em ponto flutuante. O valor "default" de Digits é 10. > Digits:=15; > cos(75.3); 90 4) Se usarmos : no final do comando, o mesmo é executado, porém o resultado não é mostrado. %; mostra o conteúdo do último output e %%; do penúltimo. > log[3](10): > %; > evalf(%%); 5) O Maple pode trabalhar com inteiros muito grandes > 253!; O comando length (%) mostra o número de dígitos de 253! > length(%); 6) Para fazer uma atribuição à variável z e posteriormente apagar o conteúdo de z: > z := 7; > z; > z := 'z'; > z ; 7) O modo texto pode ser obtido clicando em na barra de menu ou # pode ser usado para fazer comentários > # o conteúdo de z também pode ser deletado usando unassign ( 'z' ); > z :=7; > z; > unassign ( 'z' ); >z; 8) restart; reinicializa o MAPLE Exercícios: 1) Calcular o valor de : a) seno(30), seno(300) = seno ⎟ ⎠ ⎞ ⎜ ⎝ ⎛ 180 30π b) )d/cln( ba 50 35+ para a = -3,02 , b = 10, c = 10,3 e d = 7,1. Armazenar o resultado na variável v. c) secante hiperbólica de v calculado no item b. 2) O valor de a = 163e π é um número inteiro? Estimar a com 18, 29, 30, 31 e 57 dígitos. Comparar os resultados. 3) Calcular 325! pelo Maple e pela sua calculadora. Analisar os resultados. Qual o número de dígitos de 325! ? 4) O que é maior 2227 ou 2524? 5) É possível calcular 33333333333333333333 no Maple? 91 Aula de Laboratório2 utilizando Maple Objetivos : Localizar as raízes da equação f(x)=0(algébrica ou transcendente) e calculá-las utilizando comandos do Maple e o método da Bissecção. Aplicação: O fator de atrito λ para um duto retangular , segundo Maubach, é dado por: ( ) 989,0Re10log035,21 −= λ λ , onde Re é o número de Reynolds. Determinar λ para os seguintes valores de Re : 10000, 20000, 50000, 100000, 1000000. Comandos Comentários 1 f:=x->f(x); Define a função de variável x , f(x) 2 plot(f(x),x= a..b,y=c..d); Plota a função f(x) sendo: x - variável y - parâmetro na vertical(opcional) a,b,c,d :parâmetros a serem especificados 3 plot({f(x),g(x)},x= a..b,y=c..d); Plota as funções f(x) e g(x) num mesmo sistema de eixos. Para opções utilizar o Help: < ?plot 4 solve(f(x)=0, x); fsolve(f(x)=0, x); fsolve(f(x)=0, x=a ..b); Calcula as raizes da equação f(x)=0. Para opções utilizar o Help: < ?solve < ?fsolve 5 For i from io by k to in do comandos od; Comando de repetição onde: i - variável io - valor inicial in - valor final k - passo 6 if condição then comando1 else comando2 fi; Comando de teste (condicional) Exemplos 1) >f:=x->sin(x/2); > f (Pi); f (4); 2) >plot(f(x), x=-3*Pi..3*Pi); > plot(f(x), x=-15..15, y=-2..2); 3) > plot({x^2-5*x, sqrt(x+1)}, x= -2..7, y= -7..7); 4) vide exercício 1(página2) 5) > for i from 0 by 2 to 10 do > print(i^2) > od; 6) > a:=-23; > if a<0 then > print(-a); > else > print(a); > fi; 92 2) Raiz de multiplicidade ímpar > p2:=x-> 2*x^4+6.8*x^3-16.56*x^2-86.756*x-85.1690; > plot(p2(x), x=-5..5); > fsolve(p2(x)=0,x); 3) Raizes complexas > p3:=x-> x^3-2.8*x^2+8.2*x+15.6; > plot(p3(x),x=-2..5); > fsolve(p3(x)=0,x); > fsolve(p3(x),x, complex); > p4:=x-> x^4-2*x^3+11*x^2-18*x+18; > plot(p4(x),x=-3..3,y=-5..30); > fsolve(p4(x)=0,x); > fsolve(p4(x),x, complex); 4) > subs( x=2, x^2+x+1 ); 5) > diff(x^7,x); > diff(x^7,x,x); > diff(x^7,x$3); 6) Método de Newton-Raphson >f := x ->cosh(x)*cos(x)-1; > plot(f(x), x=-3*Pi..3*Pi); > plot({cos(x), 1/cosh(x)}, x=-10..10); > x[0]:=4.5; > for i from 0 to 5 do x[i+1]:=evalf(x[i]-(f(x[i])/subs(x=x[i],diff(f(x),x)))) od; > abs(x[5]-x[6]); > Digits:=20; > abs(x[5]-x[6]); 95 Aula de Laboratório4 utilizando Maple Objetivos : Resolver sistemas de equações lineares utilizando diversas opções de comandos do Maple. Analisar graficamente a solução , quando possível, no plano e no espaço. Aplicação: Equacionado o circuito resistivo mostrado na figura acima pelo método dos nós, obtemos: ⎪ ⎪ ⎪ ⎩ ⎪ ⎪ ⎪ ⎨ ⎧ = − + − + − = − + − + − + − + − + − = − + − + 0 5 2CV 1 BVCV 8 8CV 0 8 0BV 5 2BV 1 CVBV 8 8BV 5 10BV 1 AVBV 0 5 10AV 1 BVAV 8 AV Aplicando a Lei de Ohm: ⎩ ⎨ ⎧ −= −= BVCV2I BVAV1I . Determinar I1 e I2. Comandos Comentários 1 array (1..m, 1..n) Cria uma matriz cujos os elementos estão definidos.Vide help para opções 2 augment( M1, M2); Cria uma nova matriz colocando M1 à esquerda da M2 . As matrizes devem ter o mesmo número de linhas. 3 evalm(expressão matricial); Avalia uma expressão contendo matrizes. 4 gaussjord(M); Aplica o método de Gauss-Jordan 5 inverse(M); Encontra a matriz inversa de M. 6 linsolve(A,v); Calcula um vetor x que satisfaça a equação Ax=v 7 matrix(m , n, [x11,x12,...xmn]); Cria uma matriz com m linhas e n colunas com elementos x11, ...,xmn. 8 solve(eqn, var); Resolve simbolicamente equações eqn para variável var. Para opções vide help. 9 transpose(M); Determina a matriz transposta de M. 10 vector(n, [x1, x2, ..., xn]); Cria um vetor de n elementos x1, ..., xn.. 11 with(linalg); Carrega a biblioteca de álgebra linear do Maple V. 12 with(plots); Carrega o package gráfico 12 A&*S; Expressa multiplicação de matrizes (não comutativa) Exemplos : > with(plots); > with(linalg); 1) Sistema Linear Compatível Determinado (possui soluções) > solve({2*x+3*y=18,3*x+4*y=25},{y,x}); > plot({(18-2*x)/3, (25-3*x)/4}, x=0..6); > solve({3*x+2*y-5*z=8,-x+2*y+z=2,-x+2*y+3*z=4},{z,y,x}); > plot3d({(8-3*x-2*y)/(-5),2+2*x-2*y,(4+x-2*y)/3}, x=0..5,y=0..5, axes=box); 96 2) Sistema Linear Compatível Indeterminado (possui infinitas soluções) > solve({2*x+y=50,4*x+2*y=100},{y,x}); > plot({50-2*x, (100-4*x)/2}, x=0..6); > plot3d({50-2*x, (100-4*x)/2}, x=0..6,y=0..6, axes=box); > solve({3*x+y+2*z=0,-9*x-3*y-6*z=0},{z,y,x}); > plot3d({(-3*x-y)/2,2+2*x-2*y,(9*x+3*y)/(-6)}, x=-2..2,y=-2..2, axes=box); outra maneira: > solvefor[t]( x+y=1, x-y+z*t=3 ); > solvefor[x]( x+y=1, x-y+z*t=3 ); 3) Sistema Linear Incompatível (não admite solução) > solve({x+y=3, x+y=-5},{y,x}); > plot({3-x, -x-5}, x=-10..10); > plot3d({3-x, -x-5}, x=-10..10,y=-6..6, axes=box); 4) Resolver o sistema e calcular o resíduo produzido pela solução encontrada: ⎪ ⎪ ⎩ ⎪ ⎪ ⎨ ⎧ =−++ −=++− =−++ =+−+ 9.94x3x5.22x2.01x 9.34x2.53x2x1x3.0 9.214x5.83x42x5.01x4.0 7.24x3x1.02x1x2 > A := array([ [2,1,-.1,1], [.4,.5,4,-8.5], [.3,-1,1,5.2], [1,.2,2.5,-1] ]); > F := vector([2.7,21.9,-3.9,9.9]); > S := linsolve(A,F); > F1 :=evalm(A&*S); > Resíduo := evalm(F-F1); 5) Resolver o sistema utilizando o Método de Gauss-Jordan (Matrix Inversa). ⎪ ⎩ ⎪ ⎨ ⎧ =++ =++ =++ 3b3x42x32x5 2b3x22x31x 1b3x72x1x2 onde: a) b1 =16 b2 = -5 b3=11 b) b1 =25 b2 = -11 b3 = -5 c) b1 =3 b2 = 5 b3 = -5 > B := matrix(3,3, [2,1,7,1,3,2,5,3,4] ); > v1 :=vector(3,[16,-5,11] ); > v2 :=vector(3, [25,-11,-5] ); > v3 :=vector(3, [3,5,-5] ); > augment( B,v1,v2,v3 ); > gaussjord(%); 97