ResumoDoLivro

CAP.1 – ARITMÉTICA COMPUTACIONAL 1.1. Introdução O objectivo de qualquer método numérico consiste em, pela sua própria definição, fornecer  soluções numéricas a problemas matemáticos. É pois essencial compreender a noção de número, os vários tipos de números, as suas diferentes representações (sobretudo as usadas pelos computadores), as operações admissíveis, os erros cometidos e seus efeitos nos resultados. 1.1.1. A Contagem e os Números Reais Primeiro Primeiro foram os dedos. O advento da agricultura agricultura e comércio comércio impuseram métodos métodos de contagem mais sofisticados e maior domínio das operações aritméticas. Criaram-se diferentes sistemas de numeração, todos eles com um conjunto (finito e geralmente pequeno) de símbolos básicos (os algarismos) e um conjunto de regras de formação dos números. Ex. o sistema romano de numeração, as regras eram complicadas usando ora o princípio aditivo ora o subtractivo. Desta forma os algarismos podiam adquirir um valor diferente consoante a sua posição, ou seja, tinham um valor posicional. O número de símbolos básicos é designado por base desse sistema. 1.1.2. A Medição e os Números Reais Os núm número eross nat natura urais is reso resolv lver eram am o probl problem emaa da conta contagem gem,, ma mass para para me medi diçõ ções es não são são suficientes para exprimir o valor de certas grandezas físicas. Por outro lado o resultado de certas operações conduzem a resultados não inteiros. O conceito de número real permite resolver estes 2 problemas: um de natureza física e outro de natureza matemática. 1.2. Representação de Números Inteiros Um inteiro positivo N com n dígitos possui a representação decimal N = (dn dn-1 ... d1 d0)10 = dn x 10n + dn-1 x 10n-1 + ... + d1 x 101 + d0 x 100 Esta ideia pode ser generalizada a uma base b diferente de 10, em que b>=2 e inteiro N = (dn dn-1 ... d1 d0)b = dn x bn + dn-1 x bn-1 + ... + d1 x b1 + d0 x b0 Como a um número correspondem tantas representações quantas as bases, logo, infinitas representações, convém não confundir um nº com qualquer das suas possíveis representações. Um problema que se coloca é a mudança de base: Mudança da base b para a base 10 – soma ponderada Mudança da base 10 para a base b – divisões sucessivas onde se vai aproveitando os restos. Basta recordar que dividindo N por b obtemos dnbn-1 + dn-1bn-2 +...+ d1 como quociente e b0 como resto e assim sucessivamente. É usual designar por ‘bit’ (‘binary digit’) o elemento de memória básico que assume os dois estados que se associam aos dígitos 0 e 1. O número de bits disponíveis para a representação de números inteiros determina qual o maior nº inteiro representável em computador. A representação decimal desperdiça bits e é, portanto, menos económica que uma base que seja potência de 2. acresce ainda que a aritmética decimal é mais difícil de implementar em computador. 1.2.2. Inteiros Não Positivos Geralmente reserva-se o bit 0 para o sinal + e o 1 para o sinal -. Nesta representação há 2 zeros, o +0 e o –0, correspondendo portanto a 2 configurações de bits distintas. Esta duplicidade pode constituir um inconveniente, por exemplo na operação de comparação. Por isso usam-se outras técnicas que evitam este tipo de problemas e tornam mais fáceis as operações com nºs negativos. 1.3. Representação de Números Reais Na base decimal 43.82 é interpretada como 43 + .82 em que, por sua vez: 43 = 4x101 + 3x100 e .82 = 8x10-1 + 2x10-2 Generalizando: x = (dndn-1...d0.d-1d-2...d-k)10 = dnx10n + dn-1x10n-1+...+d0 + d-1x10-1+d-2x10-2 + ...+d-kx10-k Para base b é só substituir o 10 por b Os dígitos (dndn-1...d0) são a parte inteira e os (d-1d-2...d-k) a parte fraccionária

Mudança da base 10 para a base b A parte inteira é igual (divisões sucessivas); a parte fraccionária é por multiplicações sucessivas pela base para a qual queremos passar (b). Se multiplicarmos x por b verificamos que d-1 é a parte inteira do resultado e d-2b-1 + ... + d-kxb-k+1 a parte fraccionária do resultado. Se continuarmos obtemos d-2 e assim sucessivamente. É de notar que há partes fraccionárias decimais que têm representação binária finita e outros infinita, sendo que esses nºs não são pois representáveis num computador binário, que é por natureza uma máquina com capacidade finita. Mudança da Base b para a Base 10 É o mesmo método da soma ponderada pelo valor posicional do dígito na base original, a b. 1.3.1. Notação Científica de Números Reais Nas aplicações científicas há necessidade de recorrer a números muito grandes e a números muito pequenos, sendo que a maioria dos dígitos são zeros e, por outro lado, as medições que permitem obter estes nºs nem sempre garantem tantos dígitos exactos. A forma de resolver as dificuldades inerentes à representação destes nºs é a notação científica: x=+-mbt m – nº real não negativo (mantissa); b>=2 é um inteiro positivo (base); t é inteiro (expoente). Como, fixada a base b, a representação não é única, há que resolver ambiguidade -> convenção: m=0 se x=0 ; b-1 <= m < 1 se x!=0 ---> representação normalizada. Mantêm-se como ambiguidades x=0 (expoente t é arbitrário) e para números cuja mantissa tem infinitos dígitos repetindo-se periodicamente (ex: .9999... que consideramos equivalente a 0.1x101) 1.3.2. Representação em Ponto Flutuante A notação científica, tal como apresentada, não pode ser representada em computador, pois, para cobrir todos os reais, a mantissa e o expoente exigiriam um nº infinito de dígitos. Assim, é modificada no sentido de se utilizar um nº finito de p dígitos para a mantissa e um nº finito de q dígitos para o expoente, obtendo-se a representação em ponto flutuante. A notação é FP(b,p,q) , que é constituído por todos os nºs reais x da forma: x=+-mbt em que b-1 <= m <= 1-b-p |t| <= bq-1 e ainda x=0 Se o número a representar exceder o maior que se pode representar – ‘overflow’; se for menor  que o menor representável – ‘underflow’. Pi também não é representável porque exece o nº de dígitos da mantissa. 1.4. erros Na Aritmética em Ponto Flutuante 1.4.1. Erros de Representação Como acabámos de ver, a representação em ponto flutuante só permite a representação exacta de alguns nºs reais. a questão legítima é: dado um nº real x qual o nº em FP(b,p,q), que denotaremos por A(x), que o representa? Se x tiver representação exacta A(x)=x. Se não, há 2 técnicas: truncatura e arredondamento. Truncatura: Desprezam-se simplesmente os dígitos do nº real x que não cabem na mantissa, isto é, os dígitos da mantissa além dos p primeiros são desprezados. Arredondamento: O nº real x é representado pelo nº do sistema FP(b,p,q) que lhe está mais próximo em valor absoluto. A ambiguidade do arredondamento resolve-se pela técnica do arredondamento simétrico: - Se b e b/2 forem pares, arredondar de modo a que o último dígito fique ímpar; - Se b for par e b/2 for ímpar (como é o caso quando b=2 ou b=10), arredondar de modo a que o último dígito fique par; - Se b for ímpar, optar por uma das regras anteriores, já que parece não haver vantagens ou desvantagens determinantes em qualquer uma delas. Erro absoluto: E=A(x) – x = m~b t – mbt = (m~-m)bt como |m~-m| <= b -p em T e ½ b-p em A, O erro absoluto é majorado por  |E| <= b t-p em T e ½ bt-p em A Se x!=0 o erro relativo é dado por:

e=(x~-x)/x = (A(x)-x)/x - exprimindo em valores absolutos e majorando o numerador e minorando o denominador, chega à majoração para o erro relativo: |e| <= b1-p para T e ½ b1-p para A O segundo é conhecido por unidade de arredondamento do sistema – u de ponto flutuante do computador, mesmo quando este opera por T. Nestas condições é sempre válida a seguinte relação: A(x) = (1+e)x com |e| <= u (3) 1.4.2. Erros nas Operações Aritméticas Para a soma e subtracção temos: x1+- x2 = m1bt1 +- m2bt2 = (m1+-m2b-(t1-t2)) bt1 se t1>t2 (m1b-(t2-t1) +- m2) bt2 se t1<=t2 Estas operações desenvolvem-se em computador de acordo com os seguintes passos: 1. Decomposição dos operandos, isto é, separação destes números nas respectivas mantissas e expoentes; 2. Alinhamento das mantissas e dos expoentes para a soma e subtracção 3. Operação com as mantissas e/ou com os expoentes 4. Normalização da mantissas, isto é, translação à esquerda da mantissa com decremento do expoente em uma unidade por cada zero à esquerda da mantissa 5. Arredondamento da mantissa 6. Recomposição do resultado, ie, reunião da mantissa e do expoente para formar o resultado no sistema de ponto flutuante. De notar que só as operações com as mantissas introduz erro (com os expoentes só em overflow e underflow) pois são inteiros. Outra nota importante é que as operações em ponto flutuante não respeitam, em geral, as propriedades comutativa, associativa e distributiva. Soma y=x1+x2 e designando por e3 o erro relativo da soma, temos que, aplicando (3): y~= (1+e3)(1+e1)x1 + (1+e3)(1+e2)x2 (4) Esta expressão tem uma leitura importante, que é a seguinte: O resultado da soma de 2 nºs em ponto flutuante é idêntico ao que se obteria com aritmética exacta mas com operandos perturbados, isto é, com (1+e3)(1+e1)x1 em vez de x1 e... em vez de x2. As perturbações são muito pequenas, pois |e1|, |e2|, |e3| < u << 1, como sabemos. Como estas parcelas vão aparecer muitas vezes, há que estimar o seu valor: Teorema 1.4.1 – Sejam e, números tais que |ei| <= u com i=1,...,n e u um número positivo tal que nu<=0.01. Então, existe um número real θ verificando |θ| <=1 e tal que πi=1 a n (1+ei) = 1 + 1.01θnu (5) Tendo em contas o que se acabou de ver, os erros absoluto E e relativo da soma são obtidos da seguinte maneira: E=y~-y=1.01(2u)(Θ1x1 + Θ2x2) e = (y~-y)/y = 1.01(2u) ( (Θ1 (x1 / x1 + x2) + (Θ2 (x2 / x1 + x2) A majoração de e revela-se fácil: |e| <= 1.01(2u) max (|Θ1|, |Θ2|) ((x1/(x1+x2)) ((x1/(x1+x2)) + (x2/(x1+x2))) <= 1.01(2u) ou, de outra forma y~= A(x1+x2) = (x1+x2)(1+e) , com |e| <= 1.01(2u) Subtracção O tratamento da subtracção segue de perto o efectuado para a soma.

Só que a subtracção pode conduzir a erros relativos muito grandes quando os números a subtrair  são muito próximos. Este fenómeno é conhecido por cancelamento subtractivo e constitui uma fonte importante de erros nos cálculos em ponto flutuante para os quais devemos estar atentos. Multiplicação y~=A(x1x2) = [(1+e1)x1(1+e2)x2](1+e3) = (1+e3)(1+e1)(1+e2)x1x2 O que quer dizer que o resultado de y~ é idêntico ao que se obteria com operandos perturbados. E=[(1+e3)(1+e2)(1+e1)-1!x1x2 e= (1+e3)(1+e2)(1+e1)-1 e, portanto,

y~=A(x1x2) 0 x1x2(1+e)

O Teorema 1.4.1. permite estimar o erro relativo, vindo:

e = 1.01Θ(3u) com |Θ| <= 1

Divisão Repetindo a análise efectuada para a multiplicação chega-se a: (1+e1) x1 y~=A(x1/x2) = (1+e3) ---------------------(1+e2) x2 Após alguns cálculos, chegamos às expressões dos erros absoluto e relativo: E=1.01Θ(3u) x1/x2 e= 1.01Θ(3u)

com |Θ| <=1 , por conseguinte temos também para a divisão que:

y~= A(x1/x2) = x1/x2 (1+e) Somatório de Números n

y=

 xi ∑ xi i =1

Se atendermos a que o computador executa o algoritmo de somas parciais (termo a termo), vamos chegar à conclusão, após alguns cálculos, que o erro relativo pode ser facilmente majorado se os xi tiverem todos o mesmo sinal . |e| <= |Σ ηi xi/y| <= max |ηi| <= 1.01nu

com ηi=1.01Θi(n+1-i)u com |Θi|<=1

Donde se conclui que se os xi tiverem todos o mesmo sinal, o erro relativo pode ser majorado independentemente dos valores dos xi. Se não tiverem o mesmo sinal isso já não é possível, pelo cancelamento subtractivo a que já aludimos. Por outro lado, como os majorantes dos ηi decrescem com i, se for viável devemos começar por  somar os nºs mais pequenos (em valor absoluto) a fim de minimizar o produto ηixi. Embora não conduza necessariamente ao menor erro absoluto, pode contribuir para um maior precisão do resultado. Também se deve usar para os resultados parciais uma variável de dupla precisão o que conduz a um menor erro absoluto, raramente ultrapassando uma unidade de arredondamento, qualquer que seja o valor de n. Produto Interno de Vectores n

s=

∑ xiyi i =1

CAP. 2 – INTERPOLAÇÃO POLINOMIAL 2.1. Introdução 2.1.1. Generalidades Pretende-se fazer passar uma certa curva por pontos dados mas a função f tem de ter certas qualidades emergentes do contexto do problema e/ou há a necessidade de tornar o problema matematicamente bem posto. Os polinómios são excelentes candidatos devido às suas propriedades, de que iremos falar. O uso de tabelas provoca os seguintes problemas, quando o valor x não está na tabela: utilizar o valor mais próximo? E, neste caso, qual o erro cometido? As máquinas eliminaram estes problemas mas o estudo da interpolação ainda é importante devido a 2 razões: Primeira – a interpolação é um meio relativamente simples de aproximar certas funções. segunda – a interpolação constitui o fundamento de muitos outros métodos numéricos. Aos pontos x0, x1, x2, ..., xn dá-se o nome de nós de interpolação. Aos respectivos valores associados, y0, y1, y2, ..., yn, valores nodais. ||f||∞ = max |f(x)| dá-se o nome de norma (de máximo) da função f. xεΩbarra A ideia de interpolar valores pode estender-se também à interpolação de valores das derivadas, pelo que o problema fica em determinar: f (j)(xi)=yij, j=0,1,...,mi i=0,1,...,n Quando os mi forem todos iguais a zero, caímos na anterior e diz-se que a interpolação é do tipo Lagrange. Se não forem todos iguais a zero diz-se do tipo Hermite. Se apenas interpolarmos algumas das derivadas no ponto xi, dizemos que é do tipo Birkhoff.

Definição 2.1.1. Diz-se que uma função é um polinómio de grau n se puder ser escrita na forma: p(x) = anxn + an-1xn-1 + ... + a1x + a0 Aos a dá-se o nome de coeficientes, o grau é n e denota-se por deg p. O grau do polinómio nulo é, por convenção, - ∞. Se an01 o polinómio diz-se mónico. Como os computadores só fazem operações aritméticas básicas, trabalham só com funções racionais. Todas as outras funções são obtidas por meio de aproximações que envolvem polinómios. Teorema 2.1.1. (Weierstrass) Seja Ω um intervalo finito, e ٤ um número real positivo arbitrário. Então, para qualquer função f Є C(Ω), existe um polinómio p tal que ||f-p||∞ < ٤. A interpolação polinomial produz com naturalidade polinómios aproximadores. Exemplos 2.1.1. Verificar se uma dada função é um polinómio - fazer a alínea e) , dos cosenos --> polinómios de Chebyshev (de primeira espécie). 2.1.2. Algoritmo de Horner  Para calcular o valor de um polinómio num dado ponto x - p(x) -, através de um nº finito de operações aritméticas – computadores. Ãlgoritmo 2.1.1. inicialização: y=a0; w=x; para i=1 até n fazer: y=y+ai*w w=w*x fim do ciclo i p(x) = y Este algoritmo requer n somas e 2n multiplicações. Se pusermos x em evidência sucessivamente, chegamos à seguinte forma: p(x) = (((anx + an-1)x + an-2)x + ...+ a1)x + a0 Pelo que o algoritmo fica: Algoritmo 2.1.2. inicialização: y=an

para i=n-1 até 0 fazer: y=ai + y * x fim do ciclo i p(x) = y Este algoritmo é bastante melhor pois só requer n somas e n multiplicações e ainda por cima o erro de arredondamento é menor. 2.2. Formas Polinomiais A forma de representar os polinómios que vimos chama-se forma de potências simples. Há alternativas. 2.2.1. Forma de Newton p(x) = a0 + a1 (x-c) + a2 (x-c)2 + ... + an (x-c)n Ao parâmetro c dá-se o nome de centro do polinómio e a esta forma dá-se o nome de forma de potências centradas. quando c=0, a forma de potências simples é recuperada. A forma anterior  não é mais que o desenvolvimento em série de Taylor do polinómio p em torno do centro c, pelo que os coeficientes ai se podem obter por: ai = p(i)(c)/i! i=0,1,...,n A forma de Newton é obtida por generalização da expressão anterior: p(x) = a0 + a1(x-c1) + a2(x-c1)(x-c2) +...+ an(x-c1)(x-c2)...(x-cn) ci são ainda centros. O algoritmo de Horner fica então, quando aplicado a polinómios sob a forma de Newton:

Algoritmo 2.2.1 (Horner com centros) inicialização: y=an para i=n-1 até 0 fazer: y0ai + y * (x – ci+1) fim do ciclo i p(x)=y Este mesmo algoritmo de Horner permite passar de uma forma de Newton com centros ci para outra, com centros c’i Algoritmo 2.2.2. inicialização: a’n=an para i = n-1 até 0 fazer: a’i = ai + a’i+1 * (c – ci+1) fim do ciclo i p(c) = a’0 Os coeficientes aí são, como se pode provar, os coeficientes da forma de Newton. A aplicação do algoritmo uma vez permite introduzir o centro c e retirar o centro cn. A sua aplicação repetida, fazendo c=c’1, c´2, ..., c´n, sucessivamente, permite passar o polinómio na forma de potências simples para a forma de potências centradas com os centros indicados, introduzindo um novo centro de cada vez que o algoritmo é aplicado. De notar que os centros aparecem, na fórmula final, pela ordem inversa da utilizada na sua introdução, no algoritmo. Ex. 2.2.2. – Passar o polinómio p(x) = x3 – x2 + 2 à forma de Newton com centros em 1, -1 e 0. ver resolução e notas na pág.1. 2.2.2. Factorização de Polinómios A técnica de introduzir centros permite obter uma outra forma de representação bastante útil. Teorema 2.2.1. Se z1,z2,...,zk forem zeros distintos do polinómio p, então: p(x) = (x-z1) (x-z2) (x-z2) ... (x-zk) r(x) r(x) com r(x) um polinómio. Este teorema tem como corolário que o nº de zeros distintos de um polinómio de grau n é <=n. O Teorema Fundamental da Álgebra é que o nº de zeros (contando com multiplicidades) de um polinómio de grau n é n.

Teorema 2.2.2. (Unicidade do Polinómio Interpolador) Se p e q forem dois polinómios de grau <= nque assumem os mesmos valores num conjunto de nós distintos x0,x1,...,xn, então os 2 polinómios são iguais. (Demonstração...) Definição 2.2.1. Diz-se que z é um zero de multiplicidade m do polinómio p se p(z) = p’(z) = ... = p(m-1)(z) = 0 e p(m)(z) =! 0 se m=1 o zero diz-se simples, se m=2, duplo, ... Teorema 2.2.3. (generalização do teorema 2.2.1.) Se z1, z2,...,zk forem respectivamente zeros de multiplicidade m1, m2, ..., mk, então p(x) = (x-z1) m1 (x-z2)m2 ... (x-zk)mk r(x) em que r é um polinómio Teorema 2.2.4. (generalização do teorema 2.2.2. Sejam p e q dois polinómios de grau <=n e d=p-q a sua diferença. Se p e q coincidirem do seguinte modo nos nós x0,x1,...,xk d(j)(xi) = 0 j=0,1, ..., mi-1 e i=0,1,...,k com k 

∑ mi = n+1, então p=q i =0

2.3. Interpolação de Lagrange Vimos que o polinómio interpolador, quando existe, é único. O objectivo da presente secção é demonstrar a sua existência, o que pode ser feito por dois métodos. O primeiro baseia-se na formação de um sistema de equações lineares, cuja solução fornecerá os coeficientes do polinómio interpolador. Como sabemos, para que este sistema tenha solução única é necessário e suficiente que a respectiva matriz, conhecida por matriz de Vandermonde, possua um determinante diferente de zero. zero. O próxim próximoo teorema teorema vai permitir permitir provar provar sim simult ultanea aneament mentee a exi existên stência cia e uni unicid cidade ade do polinómio interpolador.

1  x0 L  1  x1 ... V(x0,x1,...,xn) =  M M O  xnL 1

v(x0,x1,...,xn) = det V (x0,x1,...,xn) =!0

( x0) n  

    M     )xnn   

( x1) n (

Teorema 2.3.1. O determinante de Vandermonde tem o valor  n

v(x0,x1,...xn) =

∏ ( xj− xi)

i , j =0   j >1

e, por conseguinte, se os nós de interpolação forem distintos, o polinómio interpolador existe e é único. A construção do polinómio por esta via tem duas desvantagens: primeira – implica resolver um sistema de equações lineares de ordem n+1 o que é de complexidade O(n3). Segunda – o sistema de equações torna-se tanto mais mal condicionado quanto maior for o grau n do polinómio, verificando-s verificando-see que este método não permite ir além de valores de n da ordem da dezena quando se trabalha em aritmética com 6 ou 7 decimais de precisão, já que os coeficientes nestas condições não têm qualquer dígito significativo. 2.3.1. Fórmula de Lagrange Uma forma mais directa baseia-se nos polinómios de Lagrange que passamos a construir. Introduzamos o polinómio nodal :

n

Wn(x) =

∏ ( x− xi) = (x-x0) (x-x1)...(x-xn) i=0

cujo grau é n+1 e que se anula em todos os nós. Por outro lado, definamos o polinómio n

lk =

∏ ( x− xi) = (x-x0)...(x-x i=0 i =!k 

) (x-xk+1) ... (x-xn)

k-1

que se anula em todos os nós menos no nó xk. O seu valor neste nó pode determinar-se por  continuidade, recorrendo à regra de Cauchy. Então o polinómio Lk(x) = lk(x)/lk(xk) anula-se em todos os nós excepto no nó xk em que toma o valor 1.

Definição 2.3.1. Os polinómios x− xi ∏ L (x) = xk− xi n

k

i=0 i =!k 

designam-se por polinómios de Lagrange relativos aos nós x0, x1, ..., xn. Recorrendo a estes polinómios, a construção do polinómio interpolador é trivial. Teorema 2.3.2. O polinómio interpolador p de grau <= n que interpola os valores y0,y1,..., yn nos nós distintos x0, x1, ..., xn é dado por: n

p(x) =

∑ L ( x) y k

k 

(Demonstração...)

k = 0

Exemplo 2.3.1. – Construir o polinómio interpolador de grau <=3 que interpola os valores seguintes: x–0 1 3 4 y - 1 -1 1 2 (ver resolução em pg.2) Apesar da sua simplicidade, a fórmula de Lagrange pode não ser a forma mais conveniente de representar o polinómio interpolador. Por 2 razões. primeira – É possível obter este polinómio realizando menos operações aritméticas. Segunda – está muito associada aos nós e estes, nas aplicações, mudam de nº ou de posição, pelo que o polinómio interpolador construído numa tentativa não é aproveitado nas tentativas seguintes. 2.3.2. Fórmula de Newton Evita inconvenientes da fórmula de Lagrange. Se considerarmos os pontos x0,x1,...,xn-1 como centros do polinómio, podemos escrevê-lo: pn = a0 + a1W0 + anWn-1 = pn-1 + anWn-1 em que W0=x-x0, ... , Wn-1 = (x-x0) ... (x (x – xn-1) xn-1) são polinómios nodais. Os coeficientes a0, a1, ... , an vão ser determinados de modo a que pn dado por esta expressão seja o polinómio interpolador nos nós x0, x1, ... , xn isto significa que devemos ter: pn8x0) = y0, pn(x1) = y1, ... , pn(xn)=yn Donde se conclui que os coeficientes do polinómio devem satisfazer: a0 = y0 ak = yk – pk-1(xk) / Wk-1(xk) Como os coeficientes só dependem dos x e dos y, se quisermos acrescentar mais pontos é só fazê-lo por indução: pn+1(x) = pn(x) + an+1Wn(x) 2.3.3. Diferenças Divididas Os coeficientes, porque dependem só dos x e dos y, podem escrever-se na forma: ak = y[x0,x1,...,xk]

estes coeficientes chamam-se habitualmente diferença dividida de ordem k e a maneira de calculá-los é estabelecida no seguinte teorema. Teorema 2.3.3. Os coeficientes ak do polinómio p de grau <= n que interpola os valores y0, y1, ..., yn nos nós distintos x0, x1, ..., xn são dados indutivamente pela expressão: ak = y[x0,x1,...,xn] = y[x1,...,xk] – y[x0,x...,xk-1] / xk – x0 Então, o polinómio pode ser escrito na forma seguinte, se rotularmos os nós de outra forma: pm,k+1(x) = (x-xm)pm+1,k(x) + (xm+k+1 – x)pm,k(x) / xm+k+1 – xm é o polinómio de grau <=k+1 que interpola nos nós distintos xm, xm+1,...,xm+k+1 os valores correspondentes dos y. Esta é a fórmula de Aitken-Neville de construção do polinómio num dado ponto x, a qual se distingue pelo facto de não precisar do cálculo específico dos coeficientes para determinar o valor  no ponto x. Da expressão do teorema 2.3.3., temos em particular que: y[x0] = y0 y[x0,x1] = y[x1] – y[x0] / x1-x0 y[x0,x1,x2] = y[x1,x2] – y[x0,x1] / x2 – x0 .... O aspecto destas expressões torna aconselhável e conveniente dispor o cálculo das diferenças divididas numa tabela (ex: para 4 nós): x y [. ] y[., .] y[. , . , .] y[. , . , . , .] x0

y0 y[x0,x1]

x1 x2 x3

y1 y2 y3

y[x0,x1,x2] y[x1,x2]

y[x1,x2,x3]

y[x0,x1,x2,x3]

y[x2,x3]

Exemplo 2.3.2. Determinar o polinómio interpolador do ex. 2.3.1. na forma de Newton (pg.3) Os valores nodais yi não tiveram até aqui qualquer ligação entre si. Mas se forem de uma funçã f, é possível estabelecer uma ligação importante entre as diferenças divididas de ordem k e a derivada da mesma ordem de f. Teorema 2.3.4. Seja f ε Ck (Ω) e x0,x1,...,xn nós distintos distintos pertencentes ao intervalo Ω. Então, existe um ponto ξ ε Ω tal que f[x0,x1,...,xn] = (1 / k! ) * f k (ξ) 2.3.4. Interpolação Inversa Se a função f possuir inversa, como é o caso quando é estritamente monótona, podemos escrever  que x=g(y) e os valores x e y podem trocar de papel --> interpolação inversa. Situação frutuosa: Supúnhamos que queríamos calcular o valor de x tal que f(x)=c. Então podemos interpolar g e depois x=g(c), ou aproximadamente p(c) Ex. 2.3.3. Determinar aproximadamente o zero da função f(x) = ln(1 + x2) – exp(-x) no intervalo [0,1]. (pg...) 2.4. Erros de Interpolação Teorema 2.4.1. Seja f ε Cn+1(Ω) e pn o polinómio de grau <=n que interpola f nos nós distintos x0,x1,...,xn contidos no intervalo Ω. Então, para qualquer ponto x pertencente a Ω existe um valor ξ ε Ω dependente de x0,x1,...,xn, de x e de f, tal que: en(x) = f(x) – pn(x) = (1 / (n+1)!) * f n+1(ξ)Wn(x)

Como o ξ é normalmente desconhecido temos de nos contentar com um majorante. assim, usando a norma do máximo é fácil de deduzir: ||en||∞ <= (1/(n+1)!) ||f n+1||∞||Wn||∞ (2.4.3.) Mas substituímos uma dificuldade por outra(s): Calcular as normas do 2ª membro. Por isso uma estimativa do erro que por vezes é útil é: ||en||∞ <= (1/4(n+1)) ||f n+1||∞hn+1 Em que h é o espaçamento máximo entre nós consecutivos. Ex: 2.4.1. Pretende-se construir uma tabela para a função f(x) 0 ln(x) no intervalo [1,2] com nós equidistantes e de modo a que o erro cometido quando se interpola linearmente nesta tabela não exceda em valor absoluto (0.5)10-5. Determinar o espaçamento dos nós. Resolução na pg.... É de notar que nos casos mais complicados, em que não conseguimos determinar ||f n+1|| ou um seu majorante aceitável, é usual recorrer a diferenças finitas (cap.3) para estimar o seu valor. Rigidez dos polinómios Uma interrogação que se põe é de saber se o erro tende para zero quando n-->infinito. a resposta é que nem sempre tal sucede. Este comportamento é por vezes referido como denunciando a rigidez dos polinómios para acompanhar a função interpolada. 2.4.2. Nós de Chebyshev A expressão 2.4.3. mostra que o erro depende de f e de W (dos pontos). A representação de W mostra que esta função possui os seus maiores valores (absolutos) na proximidade dos extremos, pelo que podemos concluir que os polinómios interpoladores devem ser usados, sempre que possível, na zona central do respectivo intervalo de interpolação. Também a extrapolação para fora do intervalo deve ser feita com cautelas. Uma questão interessante é de saber qual é a posição dos nós que tornam ||Wn|| mínimo. Veremos que é quando os nós coincidem com os zeros dos polinómios de Chebyshev. Vamos supor que o intervalo de interpolação é [-1,1]. Se não for só temos que usar a transformação de coordenadas: x = a ((1-ξ )/2) + b ((1 + ξ)/2) com ξ pertencendo a [-1,1] Teorema 2.4.2. O polinómio de Chebyshev Tn, tem os zeros localizados nos pontos (2k-1)π xk = cos ------------- , k=1,...,n 2n E os os ex extrem tremos os local ocaliizado zadoss nos nos pont pontos os x’k x’k = cos cos (kπ (kπ/n /n)) k nos quais Tn(x’k) = (-1) ----> ||Tn|| =1 Mas estes polinómios têm outra propriedades notáveis. É fácil provar que o coeficiente an de Tn é igual a 2n-1, pelo que o polinómio T’n = 2n-1Tn é mónico. Teorema 2.4.3. O polinómio T’n é, de todos os polinómios em P´n[-1,1] (todos os mónicos) o de menor norma, isto é, ||T´ ||T´n| n||| <= ||p| ||p||, |, qu qual alqu quer er qu quee seja seja p pe pert rten ence cent ntee a P´n[ P´n[-1 -1,1 ,1]] Escolhendo os zeros do polinómio de Chebyshev Tn+1 como nós de interpolação, pode sem dificuldade concluir-se que: Wn(x´k) = T´n+1(x´k) = 2 -n do nde ||Wn|| = 2-n 1 daqui resulta que en(x) = ------------------ fn+1   (ξ)Tn+1(x) n 2 (n+1)! Então a expressão do majorante do erro passa a ser: 1 ||en|| <= ----------------- ||f n+1|| em [-1,1]

2n(n+1)! (b-a)n+1 ||en|| <= ----------------- ||f n+1|| 22n+1(n+1)!

em Ω

2.5. Interpolação Com Nós Equidistantes 2.5.1. Fórmulas de Interpolação Em muitas aplicações os nós são equidistantes, pelo que há vantagem em especializar o formulário deduzido nas secções anteriores para acolher este caso. Sendo h a distância entre nós sucessivos, podemos escrever: h=(xn-x0)/n xk = x0 + kh com k=0,1,...,n É conveniente introduzir a transformação linear de coordenadas: s = s(x) = (x-x0)/h, x = x(s) = x0 + sh de que resulta f(x) = f(x0 + sh) O cálculo das diferenças divididas de uma função f pode simplificar-se se introduzirmos o conceito de diferenças conforme se explica na definição seguinte. .... 2.7 – SPLINES Mencionámos atrás, na secção 2.4., as dificuldades que podem surgir quando se empregam polinómios interpoladores de elevado grau e vimos que um remédio é escolher mais criteriosamente os nós. Outra forma de ataque é usar os splines, que são polinómios seccionalmente contínuos, uma vez que os polinómios são realmente muito atraentes em termos computacionais. Definição 2.7.1. Uma função S é um spline polinomial de grau m (m>=0) relativo aos nós a = x0 < x1 < ... < xn = b se verificar as seguintes propriedades: S coincide em cada subintervalo Ωi = [xi-1, xi), i=1,...,n com um polinómio de grau <= m S pertnece a Cm-1(Ω) O conjunto dos splines polinomiais de grau m na malha Ωi será denotado por Pm,n Empregaremos ainda a seguinte notação: hi=xi – xi-1, h = max hi Em que h se costuma dar o nome de parâmetro da malha. 2.7.1. Splines De Grau Zero é o mais simples, m=0. Coincide em cada subintervalo Ωi, com uma constante e pertence a C-1(Ω’) É óbvio que Si(x) = yi. Para construir este spline é preciso tomar uma opção relativamente à escolha dos valores dos yi’s. Se o spline interpolar uma função f ε C1(Ω) podemos escolher pontos ai Ωi e fazer Si(x)=f(ai), i=1,2,...,n. Os casos mais vulgares são tomar o extremo esquerdo do subintervalo, isto é, ai=xi-1, ou o extremo direito, isto é, ai=xi, ou o ponto médio, isto é, ai=xi-1 + xi)/2. A escolha afecta p erro de interpolação. Teorema 2.7.1. Seja f ε C1(Ω’) e S ε P0,n. Então, o erro de interpolação de f por S é majorado por  ||e|| <= ||f´|| h ou ||e|| <= ½ ||f´|| h conforme se escolham os extremos ou ponto médio. Esta expressão sugere que devemos tentar concentrar os nós nas zonas em que a primeira derivada for maior em valor absoluto, se os pudermos escolher livremente. Vê-se ainda que existe convergência quando h  0havendo vantagem em usar os pontos médios. 2.7.2. Splines de Grau 1 xi-x x – xi-1 Si(x) = yi-1 ------------ + yi yi ------------Hi hi

com xi-1 <= x <= xi

Em que os yi são os valores nodais do spline. Teorema 2.7.2. (estimativa do erro) – Seja f ε C2(Ω) e S ε P1,n. Então o erro de interpolação de f  por S é majorado por ||e|| <= 1/8 ||f’’|| h2 2.7.3. Splines Quadráticos O spline seguinte obtido pondo m=2 coincide com um polinómio de grau <= 2em cada subintervalo Ωi e S ε C1(Ω’). A construção deste spline é, contudo, menos directa. Podemos escrever que Si(x) = yi-1 + mi-1 (xi – xi-1) + mi/2 (x-xi-1) 2 (2.7.3) Em que recorremos à seguinte notação mi = S’(xi) i=0,1,...,n Mi = Si’’(x), i=1,...,n Esta forma de representação garante à partida a continuidade da derivada do spline nos nós. A expressão 2.7.3. garante que o spline quadrático assume para xi-1 o valor yi-1. Como Si’(x) = mi-1 + Mi8x – xi-1) Temos que mi = Si’ (xi) = mi-1 + Mi hi Donde se extrai Mi = (mi – mi-1) / hi, i=1,..,n (2.7.5.) Que permite obter os Mi a partir dos mi. Para que o spline interpole o valor yi em xi devemos ter  que si(xi) = yi-1 + mi-1 hi + Mi/2 hi2 = yi Combinando esta expressão com a anterior, resulta que Yi – yi-1 mi = 2 -------------- - mi-1 com i=1,...,n (2.7.6.) hi assim, se m0 for dado podemos obter através de 2.7.6. os valores nodais das derivadas m1, m2, ... e, por meio de 2.7.5. os valores dos parâmetros Mi. Contudo verifica-se que o comportamento destes splines é algo instável, razão da sua pouca utilização, havendo alternativa fácil. 2.7.4. Splines Cúbicos A facilidade de construção e a relativa estabilidade de comportamento justificam a sua maior  divulgação e popularidade. Vejamos como construí-los. Podemos escrever na forma (fica automaticamente assegurada a continuidade das 2ªs derivadas) xi-x x - xi-1 Si’’(x) = Mi-1 ------------ + Mi ----------hi hi Mi são chamados de momentos. Integrando 2 vezes: (xi – x)3 (x – xi-1)2 xi – x x – xi-1 Si(x) Si(x) = Mi-1 -------------------------- + Mi ---------------------------- + ci --------------------- + di ---------------------------- (2.7.8.) 6 hi 6 hi hi hi ci e di são constantes de integração que podem ser determinadas se impusermos as condições de interpolação: Si(xi-1) = yi-1 e Si(xi) = yi Chegando-se a: ci = yi-1 – Mi-1 hi2 / 6 e di = yi – Mi hi2 / 6 As quais se introduzidas em 2.7.8. dão: (xi – x)3 (x – xi-1)3 xi – x x – xi-1 2 2 Si(x) = Mi-1 ------------- + Mi ------------------- + (yi-1 – Mi-1 (hi /6)) ---------- + (yi – Mi (hi /6) ---------6 hi 6 hi hi hi Para concluir a construção do spline falata determinar os valores dos momentos, o que se faz impondo a condição de continuidade das primeiras derivadas dos nós, Si’(xi-) = Si+1’(xi+=, i = 1,2,...,n-1 (2.7.11) Derivando a expressão de Si(x) temos que Si’(x) = - mi-1 (xi – x) 2 / (2 hi) + Mi (x – xi-1)2 / (2hi) + (yi – yi-1)/hi – (mi – Mi-1)hi/6 E, por conseguinte Si’(xi-) = (yi – yi-1)/hi + hi Mi-1/6 + hi Mi/3

Si+1’(xi+) = (yi+1 – yi) / hi+1 – hi+1 h i+1 mi/3 – hi+1 mi+1/6 Introduzindo estas relações em 2.7.11 e agrupando termos, vem que: (hi/6) * Mi-1 + ((hi + hi+1)/3) * Mi + ((hi+1/6) * Mi+1) = (yi+1 (yi+1 – yi)/hi+1 - (yi – yi-1)/hi com I = 1,2,…,n-1 Estas expressões formam um sistema sde n-1 equações a n+1 incógnitas. Pelo que temos de impor 2 condições suplementares, que serão determinadas pela finalidade do spline Spline Completo Corresponde ao caso em que as primeiras derivadas nos nós extremos são conhecidas S1’(x0) = y0’ e S’n(xn) = y´n com y´0 e y´n dados. Então as condições suplementares são facilmente extraídas, vindo: (y1 – y0)/h1 – h1 M1/6 – h1 M0/3 = y’0 2 .7. 15. a (yn – yn-1)/hn + hn Mn-1/6 + hn Mn/3 = y’n 2 .7. 15 . b O spline assim obtido designa-se por spline completo. Spline Natural Na ausência de qualquer informação dos extremos é frequente optar pelas condiçoes S’’1(x0) = M0 = 0 S’’n((xn) = Mn = 0 Dizendo-se que neste caso o spline é natural. Estas condições fronteira reduz a precisão do spline Spline Periódico Se a função a interpolar for periódica no intervalo ómega barra, as condições suplementares devem evidentemente ser: Y0 = yn, S’(x0) = S’(xn) M0 = Mn Para assegurar a periodicidade do spline cúbico. A imposição destas condições tem como consequência nefasta a destruição da estrutura tridiagonal do sistema de equações. Continuidade Da Terceira Derivada em x1 e xn-1 Isto é, fazendo S1(3)(x1-) = S2(3)(x1*) e Sn-1(3)(xn-1-) = Sn(3)(xn-1+) Equações Resolventes Estudamos só o caso do spline completo. Para este o sistema de equações resultante das expressões 2.7.11 e 2.7.15 escrito na forma matricial é: 0    M 0    b0    2 λ 0 0 L  µ 2 λ  L    M     b   0 1 1     1    1  

M  0 0 

O

O

L

 µ n−1

L

0

   M   =  M           2 λ n −1    M n−1    bn−1         2    µ n    M n    bn   O

M

em que pusemos, por economia de notação, λ0 = 1, b0 = (6/h1) ((y1-y0)/h1 - y’0) μn = 1, bn = (6/hn) (y’n – (yn – yn-1)/hn) e, para i = 1,..., n-1 λi = hi+1 / (hi + hi+1) μi = hi / (hi + hi+1) = 1 – λi 6 yi+1 – yi yi – yi-1 bi = ------------ (---------------- - ------------- ) hi + hi+1 hi+1 hi Trata-se, como facilmente se verifica, de um sistema cuja matriz é de diagonal estritamente dominante por linhas a qual, como se demonstrará no cap.6 é sempre invertível, pelo que tem solução única qualquer que seja o 2º membro b. Acresce que é tridiagonal o que torna a obtenção dos momentos relativamente fácil, como veremos no cap.6 As boas propriedades do spline cúbico, nomeadamente a ausência de oscilações espúrias encontram justificação teórica no

Teorema 2.7.3. (Holladay) – Sejam dados os nós a 0 x0 < x1 < ... < xn = b e os valores nodais y0, y1, ..., yn. Então, de todas as funções f ε C2 (Ω’) que interpolam estes valores, o spline cúbico natural é a única função que torna mínimo o valor de J(f) =

b

∫  a

[ f ''( x))]]2 dx

J(f) representa uma espécie de ‘curvatura média’ de f no intervalo ómega barra. O enunciado diznos que o spline cúbico natural é, de todas as curvas interpoladoras com derivadas contínuas até à segunda ordem, aquela que é mais ‘direita’. Erros de Interpolação A dedução de estimativas dos erros de interpolação dos splines cúbicos está fora do âmbito. No entanto é útil saber os resultados:

Teorema 2.7.4. Seja f ε C4(Ω´) e S o spline cúbico satisfazendo qualquer das condições suplementares referidas nesta subsecção. Então, ||f – S|| <= (5/384) ||D4 f|| h4 ||D(f-S)|| <= ((sqrt3/216) + (1/24)) ||D4 f|| h3 ||D2(f-S)|| <= ((1/12) + (1/3)/(h/h´)) ||D4 f|| h2 ||D3(f-S)|| <= ½ (1 + (h/h´)2) ||D4 f|| h em que h’ é o mínimo de hi com i de 1 a n

CAPÍTULO 3 – DIFERENCIAÇÃO NUMÉRICA 3.1. Introdução Em muitas circunstâncias torna-se necessário obter valores das derivadas de uma função sem recorrer à respectiva expressão analítica por esta não ser conhecida ou por ser demasiado complicada. Por isso é conveniente ter técnicas alternativas à derivação analítica que sejam simultaneamente fáceis de usar e que permitam obter a precisão necessária. Estas técnicas são designadas por diferenciação numérica. Mas atenção, 2 funções podem convergir tanto quanto se queira e o mesmo não acontecer com as suas derivadas. Não obstante os exemplos pessimistas que se podem dar, é possível fundamentar a obtenção de derivadas substituindo a função a derivar por uma outra que de algum modo a aproxime mas cuja derivação seja mais simples. Como candidatos a funções aproximadoras surgem com toda a naturalidade os polinómios. 3.2. Derivadas De Primeira Ordem Consideremos então o problema de determinar o valor da primeira derivada de uma função f num ponto x dado, sendo conhecidos os valores nodais yi = f(xi) nos n+1 nós distintos xi = 0,1,...,n. Como referimos, a técnica mais usual consiste construir o polinómio interpolador nestes pontos e calcular analiticamente a respectiva derivada, esperando deste modo obter uma aproximação suficientemente boa. Como vimos no cap.2 en(x) = f(x) – pn(x) Com en(x) = f[x0,x1,...,xn]Wn(x) n

Com o polinómio nodal Wn(x) =

∏ ( x − x ) i

(3.2.3.)

i=0

Derivando e aplicando o teorema 2.3.4. temos uma forma de exprimir o erro da derivada esclarecedora: 1 1 (n+2) e’n(x) = ----------- f  (ξ) Wn(x) + -------------- f (n+1) (η) w’n(x) ( 3 . 2 .5 . ) (n+2)! (n+1)! onde ξ e η pertencem ao intervalo ómega barra. Em duas situações especiais mas com interesse prático, esta expressão pode ser simplificada O ponto x coincide com um dos nós Tendo em atenção 3.2.3., o 1º termo de 3.2.5. anula-se vindo: n

∏

e’n(x) = (1/(n+1)! )* f (n+1) (η) W’n(xi) = (1/(n+1)!) f (n+1) (η)   j = 0 ( xi − x j ) =j!

i

O ponto x é um zero de W´n Em certos casos o ponto em que se pretende calcular a derivada está disposto simetricamente em relação aos nós, o que implica imediatamente que estes sejam em nº par e que, portanto, n seja ímpar. Se designarmos agora por x barra o ponto em que se quer calcular a derivada, devido à simetria temos que W’n(x) = 0 e então o 2º membro de 3.2.5 anula-se e obtemos: ( n−1) / 2

(n+2)

e´(x) = (1/(n+2)!) f 

(ξ)

∏ (−( x&&&− x ) ) 2

i

i =0

Vejamos agora alguns casos particulares de maior aplicação 3.2.1. Diferenças Finitas De Primeira Ordem n=1 e x=x0 p1(x) = f(x0) + f[x0,x1](x-x0) p’1(x) = f[x0,x1] e, portanto, f’(x0) ≈ p’1(x0) = f[x0,x1] = f(x1) – f(x0) / x1 – x0 pondo h = x1 – x0 e usando Dhf(x) para indicar a derivada aproximada de f no ponto x dhf(x) = f[x + h, x] = f(x+h) – f(x) / h Esta expressão é a diferença diferença finita progressiva de primeira primeira ordem e h é o passo.

Como x é um nó podemos aplicar 3.2.6. 3.2.6. para obter o erro : e’1(x) = ½ hf’’(ξ) com ξ ε [x0,x1] Daqui decorre que se a 2ª derivada de f é contínua em [x0,x1], o erro tende para zero com a primeira potência de h N=1 e x=x1 Da mesma maneira obtém-se a diferença regressiva de primeira ordem e o erro: Dhf(x) = f[x,x-h] = f(x) – f(x-h) f(x-h) / h e’1(x) = ½ hf’’(ξ) n=1 e x=(x0+x1)/2 Chega-se sem dificuldade a Dhf(x) = f[x+h,x-h) = f(x+h) – f(x-h) / 2h Que é a diferença finita central de primeira ordem O erro é e’1(x) = -1/6 h2f (3)(ξ) Para funções em que a 3º derivada é contínua no intervalo, o erro tende para zero com h2, o que representa uma melhoria. Exemplo 3.2.1 (fazer) – Calcular a primeira derivada da função f(x) = exp(sin x) no ponto x=0,5 utilizando as diferenças finitas de primeira ordem com h=0,01 3.2.2. Diferenças Finitas de 2ª Ordem Todas as fórmulas anteriores foram baseadas em interpolação linear. Vejamos agora se os polinómios forem do 2º grau, isto é, n=2 p2(x) = f(x0) + f[x0,x1](x-x0) + f[x0,x1,x2](x-x0)(x-x1) p’2(x) = f[x0,x1] + f[x0,x1,x2](x-x0 + x-x1) X=x0 Dhf(x0) = f[x0,x1] + f[x0,x1,x2](x0-x1) e e’2(x0) = 1/6 f (3)(ξ) (x0-x1)(x0-x2) Se os nós forem equidistantes: -3f8x) + 4f(x+h) – f(x+2h) Dhf(x) = ------------------------------------------2h 2 (3) E e’2(x) = 1/3 h f  (ξ)

( 3 . 2 .1 4 )

( 3 . 2. 1 6 )

3.2.14 e 3.2.16 são conhecidas por diferenças finitas progressivas de 2ª ordem. Uma observação pertinente é que apesar de usarem 3 pontos de interpolação apresentam uma estimativa de erro pior que a central de 1ª ordem – simetria é benéfica. x=x2 chegamos a Dhf8x2) = f[x0,x1] + f[x0,x1,x2] (2x2 – x0 – x1) e’2(x) = 1/6 f (3)(ξ)(x2-x1) (x2-x0) 3f(x) – 4f(x-h) + f(x-2h) e a Dhf(x) = --------------------------------------------------------------------2h 2 (3) e’2(x) = 1/3 h f  (ξ) x=x1 Dhf(x1) = f[x0,x1] + f[x0,x1,x2] (x1-x0) e’2(x) = 1/6 f (3)(ξ)(x1-x0)(x1-x2) Se os nós forem equidistantes, dará: Dhf(x) = f(x+h) – f(x-h) / 2h Que é, nada mais nada menos, a diferença central de 1ª ordem 3.3. Derivadas de Segunda Ordem As técnicas anteriores generalizam-se facilmente ao cálculo de derivadas de ordem superior. Agora, é claro que o primeiro polinómio a considerar é de grau 2, isto é, n=2 interpolando nos nós x0,x1 e x2:

p2(x) = f(x0) + f[x0,x1](x-x0) + f[x0,x1,x2](x-x0)(x-x1) p´´(x) = 2f[x0,x1,x2] e, por conseguinte, f’’(x) ≈ p’’(x) = 2f[x0,x1,x2] Se os nós forem equidistantes: D2hf(x) = f(x-h) – 2f(x) + f(x+h) / h2 Derivando 2 vezes: e’’2(x) = f[x0,x1,x2,x,x,x]W2(x) + 2f[x0,x1,x2,x,x]W’2(x) 2f[x0,x1,x2, x,x]W’2(x) + f[x0,x1,x2]w’’2(x) Se x coincidir com algum dos nós resulta que W2(x)=0 e, substituindo as diferenças divididas por  derivadas temos e’’2(x) = 1/12 f (4)(ξ)W’2(x) + 1/6 f (3)(η)W’’2(x) No caso de os nós serem equidistantes: e’’2(x0) = 1/6 h2f (4)(ξ) – h f (3)(η) e’’2(x2) = 1/6 h2f (4)(ξ) – h f (3)(η) x=x1 - h anula-se e: e’’2(x1) = -1/12 h2f (4)(η) Onde mais uma vez se mostra a vantagem das diferenças finitas centrais. A tabela da página 105 compila todas as fórmulas deduzidas e mais algumas. Fórmula geral n

(k)

k

f  (xj) ≈ k!/n! * 1/h

∑ a f (x ) i

i

i =0

Exemplo 3.3.1. (fazer) – Uma partícula move-se sobre o eixo dos x, tendo-se registado as seguints posições ao longo do tempo t: t 0 ,0 0, 2 0 ,4 0, 6 0 ,8 1 ,0 x 0 ,0 0,1 987 0 ,3 8 9 4 0, 564 6 0 , 7 17 4 0 , 8 4 1 5 Pretende-se determinar a aceleração nos instantes t=0,0 e t= 0,6 3.4. Derivação Com Splines Com splines é possível obter fórmulas de diferenciação numérica de elevada precisão. Mas isso pode ser computacionalmente caro, pelo que há, em cada caso, de ponderar o custo vs. Precisão pretendida. (xi – x)2 (x – xi-1)2 yi – yi-1 S’i(x) = -Mi-1 ------------- + Mi --------------- + ------------- + (Mi-1 – Mi) hi/6 2hi 2hi hi x1 – x x – xi-1 S’’i(x) = Mi-1 ----------- + Mi -------------hi hi 3.5. Influência Dos Erros de Arredondamento As expressões de erro deduzidas mostra, que, se as derivadas apropriadas da função f forem limitadas, a derivada calculada pelas fórmulas de diferenças finitas converge para o valor exacto com uma certa potência de h. No entanto, se tivermos em conta os erros de arredondamento, esta situação altera-se radicalmente. Se utilizarmos um exemplo (pg.107), vê-se que o erro diminui numa primeira fase com h mas depois aumenta. Se calcularmos o erro total (no exemplo para as centrais) vamos obter: e. (h+h) – e.(x-h) h2 E = ----------------------- - ---- f (3)(ξ) com e. Vindo Vindo de f(x) = f.(x) + e.(x) e f. Valor aproximado de f  2h 6 Vê-se que é composto por duas partes distintas, uma proveniente da fórmula das diferenças finitas e proporcional a h2 e outra resultante dos erros de arredondamento e proporcional a 1/h.

Daqui pode inferir-se que deve haver um valor de h para o qual o E é mínimo. Infelizmente não são conhecidos nem a função e nem o valor de ξ. No entanto é possível estimar um valor  aproximado, razoável. Se admitirmos as hipóteses: ||f (3)|| <= M3 e |e.(x+h) – e.(x-h) <= 2ε Em que M3 designa um majorante das terceiras derivadas e ε um parâmetro que depende da forma como a função f é calculada. Ε será da ordem de grandeza da unidade de arredondamento do computador: |E| <= ε/h + M3/6 * h2 O 2º membro da expressão é mínimo quando h = (3ε / M3)1/3 3.6. Extrapolação de Richardson Vimos que os erros de arredondamento limitam seriamente a precisão que é possível obter com as fórmulas de derivação numérica. Vamos ver que existe um processo de minorar, mas não eliminar, o efeito da aritmética de precisão finita. Consideremos a fórmula de diferenças finitas progressivas de primeira ordem, para concretizar. O erro pode ser expresso por: ∞

2

e(x) = f´(x) – Dhf(x) = f’(x) – (f(x+h) – f(x))/h = -1/2 f’’(x)h – 1/3! f’’’(x)h - ... =

∑c h

i i

i =1

ou seja, o erro é dado por uma série de potências de passo h em que os coeficientes ci envolvem as derivadas f (i), sendo, por isso difícil de calcular em geral. Nestas condições ∞

Dhf(x) = f’(x) +

∑c h

i

i

(3.6.2.)

i =1

O processo de extrapolação de Richadson consiste em considerar a derivada dada, não por  3.6.2., mas por uma série truncada a partir de um certo termo, k 

Dhf(x) ≈ f’(x) +

∑c h

i

i

i =1

e ajustar um polinómio em h de grau <= k aos valores obtidos para uma sucessão decrescente {hi}ki=0 . O valor deste polinómio no ponto h = 0 fornece naturalmente uma aproximação, em princípio melhor, para f´(x). O facto de h=0 não pertencer ao intervalo dá o nome. Designemos por Dm,k o valor extrapolado usando os seguintes nós de interpolação hm, hm+1,..., hm+k A fórmula de Aitken-Neville permite obter este valor sem necessidade de construir explicitamente o polinómio interpolador, vindo: hm+k+1 Dm,k - hm Dm+1,k Dm,k+1 = --------------------------------------------------------------------Hm+k+1 - hm Uma sucessão frequentemente usada é hi = h0/2i, i=0,1,2,..., ou seja, procedendo por bissecções sucessivas de um intervalo original h0. Exemplo 3.6.1. (fazer) – aplicar o processo de extrapolação de Richardson ao cálculo de f’(0) com f(x)=exp(x) com diferenças finitas centrais e bissecção dos passos.

CAPÍTULO 4 – INTEGRAÇÃO NUMÉRICA 4.1. Introdução Designaremos de um modo geral por integração numérica o processo de obter valores aproximados para o integral de uma função f no intervalo Ω´= [a,b] da recta real, ou seja, b

∫ 

aproximações para I(f;Ω’) =

f ( x) dx dx

a

A necessidade de recorrer a métodos aproximados provém das seguintes situações: - A expressão analítica de f não é conhecida – é o que acontece quando esta função é dada por  tabelas ou obtida por medições de grandezas. - A expressão analítica de f é dada mas a primitiva desta função: - não é conhecida e, portanto, a forma usual de determinação do integral não é viável, ou é ‘demasiado’ complicada e, portanto, a forma usual de determinação do integral não se revela económica. Como na diferenciação, a chave do problema é aproximar a função por outra cujo integral seja fáil de calcular. Mais uma vez vamos pelos polinómios interpoladores. É razoável esperar que o valor de Ih(f) = I(pn) seja, sob certas condições, aproximado de I(f). O erro é: Eh(f) = I(f) – ih(f) = I(f) – I(pn) = I(f-pn) Como se vê, o erro depende da maior ou menor aproximação do polinómio pn a f e adiante apresentaremos estimativas desta importante grandeza.

Mudança do Intervalo de Integração Muitas vezes é necessário aplicar uma mudança de variável no integral; para transformar o intervalo [a,b] de integração num dos intervalos normalizados ou de referência usuais [-1,1] ou [0,1]. b

β 

b

α 

∫   f  ( x)dx = ∫   f  (T (ζ  )) J (ζ  ) d ζ  

Então I =

Em que J(ξ ) = T’(ξ ) é o jacobiano da transformação. Registamos, para uso futuro, as fórmulas de transformação aplicáveis neste caso: Se [α ,β ] = [-1,1]

T(ξ  J(ξ  T(ξ  J(ξ 

Se [α ,β ] = [0,1]

) = a(1 - ξ ) / 2 + b (1+ξ ) / 2 ) = (b-a)/2 ) = a (1-ξ ) + bξ ) = b-a

4.2. REGRAS BÁSICAS Seja, como habitualmente pn o polinómio de grau ≤ n que interpola a função f nos nós x0,x1,...,xn: n

pn(x) =

 xi ) Li  Li ( x) ∑  f  ( xi i =0

Sendo assim, é fácil ver que: b

Ih (f) = I(pn) =

b

n

∫ p ( x) dx = ∑ a f ( x ) = ∫ L ( x) dxdx n

a

i

i =0

i

i

a

b

Pondo Ai =

∫ L ( x) dxdx i

podemos escrever que:

a

n

Ih(f) =

∑ A f (x ) i

i

(4.2.2.)

i =0

Esta expressão costuma designar-se por regra de integração ou fórmula de quadratura , e os Ai por coeficientes ou pesos dessa regra. Consoante o valor de n e a localização dos nós no intervalo ómega barra, assim se obtêm diferentes regras de integração.

Definição 4.2.1. – Uma regra de integração diz-se de grau ( de exactidão) n se integrar  exactamente todos os polinómios de grau  ≤  n e existir pelo menos um polinómio de grau n+1 que não é por ela integrado exactamente.

Uma conclusão é que o grau de exactidão da fórmula 4.2.2. é ≥ n.

Teorema 4.2.1. – O grau de exactidão de uma regra de integração I h da forma (4.2.2.) que utilize os n+1 nós distintos x0,x1,...,xn não pode exceder 2n+1.

Uma questão que tem interesse é a de saber se este valor limite para o grau de exactidão pode ser realizado por alguma regra. A resposta é sim.

4.2.1. DEDUÇÃO DAS FÓRMULAS Correspondem a diferentes escolhas dos polinómios: 1 – REGRAS DO RECTÂNGULO caso mais simples – polinómio de grau 0, que interpola a função f no ponto x0. p0(x) = f(x0) b

e, portanto,

Ih(f) = I(p0) =

∫  f ( x ) dx = (b-a) f(x ) 0

0

a

Se fizermos coincidir x0 com o extremo esquerdo do intervalo, isto é, x0=a, obtemos a regra do rectângulo à esquerda:

Ih(f) = (b-a) f(a) Se tomarmos x0=b obtemos a regra do rectângulo à direita Ih(f) = (b-a)f(b) É fácil verificar que ambas possuem grau zero. 2 – REGRA DO PONTO MÉDIO Uma outra escolha natural será fazer x0 = (a+b)/2, obtendo-se a regra do ponto médio : Ih(f) = (b-a) f((a+b)/2) Esta regra é de grau 1. 3 – REGRA DO TRAPÉZIO Passamos agora aos polinómios interpoladores de grau n=1 cuja forma de Newton é p1(x) = f(x0) + f[x0,x1](x-x0) Teremos que… desenvolvendo obtemos: Escolhendo x0=a e x1=b, obtemos a regra do trapézio : Ih(f) = (b-a)/2 * [f(a) + f(b)] O grau é 1. 4 - REGRA DE SIMPSON Passamos agora a polinómios de grau 2. O polinómio interpolador da função f nos pontos x0,x1,x2 é: P2(x) = f(x0) + f[x0,x1](x-x0) + f[x0,x1,x2](x-x0)(x-x1) f [x0,x1,x2](x-x0)(x-x1) Tomando x0=a , x1=(a+b)/2 e x2=b e desenvolvendo: Ih(f) = I(p2) = (b-a)/6 * [f(a) + 4f((a+b)/2) + f(b)] Que é a célebre regra de integração de Simpson . O grau é pelo menos 2 sendo na realidade 3. 5 – REGRAS DE NEWTON-COTES Todas as regras apresentadas recorrem a polinómios interpoladores da função integranda em nós equidistantes no intervalo de integração. Estas regras podem ser deduzidas de forma sistemática e constituem a família das regras de Newton-Cotes: Tabela 4.2.1. – Fórmulas de Newton-Cotes n

Ih(f) = (b-a)/2 * (

∑ a f (x ) ) i

i

i =0

Xi = a + ih, I=0,1,…,n

h=(b-a)/n

Nota: Os coeficientes ai são simétricos, isto é, an-i=ai N d a0 a1 1 2 1 2 6 1 4 3 8 1 3 4 90 7 32 5 2 88 19 75 6 8 40 41 216 7 1 7280 75 1 35 77 8 2 8350 98 9 58 88 A partir de n=8 aparecem pesos negativos o que é arredondamento – promove cancelamento subtractivo.

a2

a3

a4

12 50 27 272 1 323 298 9 - 928 1 04 9 6 - 454 0 nocivo do ponto de vista dos erros de

4.2.2. ERROS DE INTEGRAÇÃO Para poder escolher a regra de integração a utilizar. Do que vimos em 2.4.1. podemos deduzir que o erro de integração é dado por: Eh(f) = I(f) – Ih(f) =

b

∫  f [ x , x , ...x , x]W ( x) dx 0

a

1

n

n

Há 2 situ situaç ações ões típ típic icas as e com com inter interes esse se práti prático co em que esta esta expre express ssão ão pod podee sofre sofrerr um umaa simplificação notável. A primeira é quando o polinómio Wn não muda de sinal no intervalo ómega barra. Invocando o teorema do valor médio para integrais b

∫ 

Eh(f) = f[x0, x1,…, xn,η ] Wn ( x)dx a

com η ∈ ómega barra. Admitindo que f ∈ Cn+1 (Ω ), por força do teorema 2.4.1., vem que:

1

b

Eh(f) = ------------- f  ( ) ∫ W (n+)

n

( x)dx

a

(n+1)! Uma conclusão imediata é que o grau de exactidão das respectivas regras é igual a n. A segunda situação é quando b

∫ W ( x)dx = 0 n

a

tendo em atenção a definição e a propriedade de simetria das diferenças divididas ... se for possível escolher xn+1 de modo a que Wn+1 possua um único sinal em ómega barra fazemos o mesmo que antes e se f ∈Cn+2(Ω ), então podemos concluir que:

1

b

Eh(f) = ---------------- f (n+2)( ) ∫ W + ( x)dx n 1

a

(n+2)! O grau de exactidão é n+1 Apliquemos estas estimativas de erro às regras de integração mais usuais. Para as regras do rectângulo temos que n=0 com W0(x) = x-a (esquerda); W0(x) = x-b (direita). b

Eh(f) = f’(ξ )

∫ 

( x − a) dx =

a

direita.

1 2

f '(ξ )( b − a) 2 , o mesmo resultado dá para a regra do rectângulo à

Para a regra do ponto médio vem que W0(x) = x-x0 em que x0 = (a+b)/2, que muda de sinal em ómega barra. Contudo o integral (segunda hipótese/variante) é zero. Escolhendo x1=x0 resulta que W1(x) = (x-x0)2 que é de sinal único. Aplicando 4.2.13 vem: Eh(f) = 1/24 f’’(ξ )(b-a)3 Tem pois grau de exactidão igual a 1. Passemos à regra do trapézio. Como netse caso W1(x) = (x-a)(x-b) tem um único sinal: Eh(f) = - 1/12 f’’(ξ )(b-a)3 também é de grau 1.Uma comparação não favorece a regra do trapézio, já que esta necessita de 2 valores da integranda, contra apenas 1 da do ponto médio para uma estimativa estimativa de erro que até é pior. No entanto, nas compostas, a do trapézio pode ser melhor. A de Simpson, temos que W2 muda de sinal, mas ... igual a zero; fazendo x3=x1=a+b/2, esulta um W3 que é de sinal único. Calculando, chegamos a Eh(f) = - 1/2880 f (4)(ξ )(b-a)5 … 4.5. REGRAS COMPOSTAS A busca de maior precisão passa pela utilização de fórmulas de grau maior. Mas isso nem sempre é assim, pois a função pode não ter a regularidade necessária. Uma alternativa são as regras compostas. A ideia é: observando as expressões do erro das várias fórmulas, verificamos que eme todas elas este depende de uma certa potência do comprimento b-a. Então, uma redução do comprimento do intervalo tenderá a reduzir o erro tanto mais quanto maior for o expoene de b-a. Então subdivide-se o intervalo em N subintervalos. Depois calcula-se o integral numérico de cada subintervalo e faz-se o seu somatório. vejamos para subintervalos todos do mesmo tamanho h= (b-a)/N Regras do rectângulo compostas:  N 

Ih(f;Ω ) = h

∑  f (a − 1) i

i =1

O erro é a soma dos erros cometidos em cada subintervalo, que vai dar:  N 

Eh(f;Ω ) = h /2 2

∑  f  '(ξ  )

mas adimitindo que f ∈C1(Ω ) e invocando o teorema do valor médio para

i

i =1

somatórios do valor de uma função:  N 

∑  f  '(ξ  ) = N f’(ξ ) = f’(ξ ) b-a/h, i

ξ∈Ω

i =1

pelo que se obtém: Eh = (b-a)/2 f’( )h Um raciocínio idêntico leva à conclusão que para a direita é tudo igual mas simétrico. Vê-se que as regras levam o erro a tender para zero com h. Perde-se uma unidade do expoente de h, ao passar das simples para as compostas(todas). Regra do Ponto Médio Composta  N 

Ih(f) = h

∑ f( a i =1

i −1

+ h / 2)

Eh(f) = (b-a)/24 f’’(ξ )h2 convergência quadrática em h para funções pertencnetes a C2 de ómega barra. Regra do Trapézio Composta  N 

Ih(f) =

∑ i =1

)] = h[1/2 f(a) + h / 2[ f (ai −1 ) + f (ai )]

 N −1

∑ f ( a ) + 1/ 2 f ( b)])] i

i =1

Eh(f) = - (b-a)/12 (b-a)/12 f’’(ξ )h2 Já não é tão desvantajosa em relação à do ponto médio composta, como acontecia com as simples. Regra de Simpson Composta  N −1

Ih(f) = h/6 [f(a)+f(b) + 2

∑ i =1

Eh(f) = b-a/2880 f  (ξ )h (4)

4

f ( ai ) + 4

N 

∑ f( a

i −1

i 01

+ h/ 2)]

Regra do Trapézio Corrigida Composta  N −1

Ih(f) = h[1/2 f(a) +

∑ f ( a ) + 1/ 2 f8 b)] + h /12 [f’(a) – f’(b)] 2

i

i 01

A desvantagem do cálculo de derivadas da função f fica agora reduzida aos extremos do intervalo e a correcção introduzida permite ter xpressão de erro mais favorável. Eh(f) = b-a/720 f (4)(ξ )h4 Uma aplicação sagaz às funções f periódicas, conduz a:  N 

Ih(f) = h

∑  f (a ) i

i =1

CAP. 5 EQUAÇÕES (algébricas) NÃO-LINEARES NÃO-LINEARES 5.1. INTRODUÇÃO 5.1.1. Raízes e Erros O problema que iremos resolver pode resumir-se deste modo: determinar os valores de z que tornam nulo o valor da função f, ou seja, resolver a equação f(z)=0. Centrar-nos-emos nas funções reais de variável real. exigiremos ainda que f possua regularidade suficiente, variando esta consoante os métodos. Exemplo 5.1.1. As funções seguintes, com domínios e regularidades diversos, são representativas do tipo de funções não-lineares cujos zeros podem ser calculados pelos métodos (aproximados) a desenvolver neste capítulo. f(x) = sin x – x2 f ε C∞( R ) f(x) = x4 – x2 + 2x1/2 f ε C (R+) f(x) = exp x - |x|3/2 - 2 f ε C1 (R)

Definição 5.1.1. Se f(z)=0, então diz-se que z é uma raiz da equação f(x)=0 ou que z é um zero da função f. Os zeros podem ser simples (a função cruza o eixo dos xx e muda de sinal), duplos (a função toca o eixo dos xx mas não muda de sinal) ou triplos (a função cruza o esixo dos xx e muda de sinal, mas fá-lo de forma a que o zero seja um seu ponto de tangência com o eixo dos xx. Definção 5.1.2. A multiplicidade de um zero z da função f é o supremo dos valores k tais que:

lim −>∞ k 

| ek  + 1 | | ek  | p

=c

se m=1, o zero diz-se simples, se m=2 diz-se duplo,... A multiplicidade não é, contudo, necessariamente, um nº inteiro. Exemplo 5.1.2. Determinar a multiplicidade de zeros (fazer): z=1; f(x) = x2-1 simples 2 2 z=1; f(x) = x – 2x + 1 = (x – 1) duplo 1/2 z=0; f(x) = x m= ½ 1/2 z=1; f(x) = x –1 simples Agora podemos compreender melhor a afirmação coloquial de que na vizinhança de um zero z de multiplicidade m a função f se ‘comporta’ como (x-z)m . Por outro lado, se f tiver uma regularidade suficiente, é possível relacionar a multiplicidade dos zeros com as derivadas desta função.

Teorema 5.1.1. Se z for um zero da função f, e se f for m vezes continuamente diferenciável no ponto z, então a multiplicidade de z é m sse f(z) = f’(z) = ... = f (m-1)(z) = 0 e f (m)(z) =! 0 5.1.2. Iterações e Ordem de Convergência Os métodos que vamos recorrer assumem, em geral, o carácter iterativo. Concretamente, estes métodos partem do conhecimento de s valores aproximados x0, x1, x2, ... , xs-1 da raiz z e com estes constróem uma nova aproximação xs desejavelmente melhor. xk = gk(xk-1,...,xk-s), com k = s, s+1, ... (5.1.1.) em que gk é uma função apropriada específica do método em questão e conhecida por função de iteração. Se esta função gk não depender de k, o método iterativo diz-se estacionário. caso contrário é não-estacionário . O esquema iterativo 5.1.1. gera, quando aplicado repetidamente, uma sucessão de valores xk aos quais correspondem correspondem os erros ek = z - xk Obviamente, é de interesse que os métodos iterativos sejam convergentes: lim xk = z lim ek = 0 k∞ k∞ Por outro lado é preciso ter em atenção a rapidez  de convergência.

Definição 5.1.3. Sejam xk e x’k duas sucessões que convergem para o mesmo limite z. Diz-se que x’k converge mais rapidamente que xk se:  x 'k  − z 

=0 lim −>∞  x − z  k 

k 

Definição 5.1.4 a) Se o erro de um método iterativo na obtenção do zero z da função f satisfizer a condição |ek+1| m <= ------------ <= M para qualquer k >=N |ek|p com 0< m <= M < ∞ , p>=1 e N um nº natural, diz.-se que p é a ordem de convergência do método relativamente ao zero da função f. b) Se, além disso, existir uma constante c>0 tal que:

lim −>∞ k 

| ek  + 1 | | ek  | p

= c diz-se que c é a constante de erro asimptótico

Se p=1, a convergência diz-se de primeira ordem ou linear, se p>1, supralinear, se p=2, de segunda ordem ou quadrática,... p e c podem depender da função e do zero também, isto é, um mesmo método pode exibir  diferentes ordens de convergência e constantes de erro em diferentes zeros duma mesma função.

Teorema 5.1.2. Seja ek uma sucessão que satisfaz a condição a) da definição 5.1.4. Se p = 1 e M<1, ek converge para zero qualquer que seja o valor de eN. Se p>1, ek converge para zero para qualquer valor de M desde que o valor eN seja suficientemente pequeno. A necessidade de ter ou não estimativas iniciais ‘suficientemente próximas’ do zero pretendido leva-nos a distinguir entre métodos de convergência local (exigência requerida) e métodos de convergência global . 5.2. MÉTODO DA BISSECÇÃO 5.2.1. Descrição do Método O método da bissecção consiste em construir subintervalos Ik = [ak,bk] ⊂ I = [a,b] por divisões sucessivas a meio e relativamente aos quais também se verifique que f(ak) e f(bk) tenham sinais diferentes. Consideremos as seguintes 3 situações. Se f(xk+1) = 0, então xk+1 é um zero e o algoritmo termina aqui. Se f(ak) e f(xk+1) tiverem sinais diferentes, faça-se ak+1 = ak e bk+1=xk+1 e ... ao contrário. 5.2.2. Erros e Convergência Teorema 5.2.1. Seja f ε C[a,b] e tal que f(a) e f(b) tenham sinais diferentes. Então, a sucessão {xk} determinada pelo método da bissecção converge para um zero de f neste intervalo, e o erro satisfaz a relação |a-b| |ek| <= ---------( 5 . 2 .2 . ) k 2 Da expressão 5.2.2. deduz-se imediatamente que |ek+1| 0 <= ------------ <= ½ |ek| Exemplo 5.2.1. (fazer) Do exemplo anterior, podemos concluir que o método da bissecção é global; A convergência poderá ser bastante lenta, pelo que, por vezes, é usado como fase preparatória para outros métodos. 5.3. MÉTODO DA FALSA POSIÇÃO

5.3.1. Descrição do Método Procede-se do mesmo modo que no método da bissecção excepto que o ponto xk se determina não como um ponto médio do intervalo, mas como a intersecção da secante que passa pelos pontos (ak, f(ak)) e (bk, f(bk)). A equação desta secante pode escrever-se nas seguintes formas alternativas: y = f(ak) + f[ak, bk](x – ak) y = f(bk) + f[ak,bk](x-bk) pelo que xk+1 é dado por uma das expressões seguintes: f(ak) xk+1 = ak - -----------f[ak,bk] f(bk) xk+1 = bk - --------------f[ak,bk] xk+1 = ak f(bk) - bk f(ak) / f(bk) – f(ak) (5.3.5.) O denominador da expressão 5.3.5. é constituído por 2 parcelas com o mesmo sinal, pelo que não é de recear o aparecimento do cancelamento subtractivo. Método da Falsa Posição Modificado É lenta quando um dos extremos, ak ou bk, se imobiliza. Então pode usar-se a seguinte modificação: empregar os pontos (ak, f(ak)) e (bk, f(bk)/2) para traçar a secante em vez dos pontos anteriores -- É o algoritmo de Illinois 5.3.2. Erros e Convergência Teorema 5.3.1. Seja f uma função convexa ou côncava no intervalo I 0 [a,b] com f(a) e f(b) de sinais diferentes. Então, o método da falsa posição converge linearmente para o zero de f neste intervalo. 5.4. MÉTODO DA SECANTE 5.4.1. Descrição do Método Consiste em, partindo de 2 iterações quaisquer (xk-1, f(xk-1)) e (xk, f(xk)), obter o valor seguinte xk+1 como a intersecção da secante que passa pelos referidos pontos com o eixo dos xx. assim, este método apenas difere do da falsa posição por não se exigir que os valores da função nos extremos dos subintervalos tenham sinais diferentes. Isto pode fazer com que este método deixe de convergir em certas situações. Todavia, quando converge, fá-lo mais rapidamente (ordem de convergência superior a 1). A expressão que permite, obter xk+1 é: yk xk-1 – yk-1 xk xk+1 = --------------------------------------------- em que y k = f(xk) yk – yk-1 Uma expressão equivalente mas menos sensível a erros de arredondamento é: yk xk+1 = xk - ----------------f[xk-1, xk] Esta pode ser escrita de uma forma equivalente mas mais elucidativa: ~ yk xk+1 = xk + hk com hk = - --------------f[xk-1, xk] 5.4.2. Erros e Convergência Teorema 5.4.1. Se todas as iterações estiverem contidas numa vizinhança (a,b) suficientemente pequena do zero z da função f ε C2[a,b], então o método da secante é convergente, e o erro satisfaz a relação |ek+1| <= M|ek| |ek-1| com M = M2 / 2m1, 0 < m1 <= |f’(ξ )| < M1 e |f’’(ξ ) <= M2 para todo o ξ ε [a,b], e a ordem de convergência é p = (1+ V5)/2 = aprox. a 1,618

5.5. MÉTODO DE NEWTON 5.5.1. Descrição do Método A curva é aproximada pela sua tangente, e a intersecção desta com o eixo dos xx é tomada como o novo valor da aproximação ao zero de f. A equação da tangente à curva y=f(x) que passa pelo ponto xk é: y = f(xk) + f’(xk)(x-xk) e, portanto, a sua intersecção com o eixo dos xx ocorre na posição f(xk) xk+1 = xk - --------------f’(xk) Este método implica que tenhamos de programar, além da função f, a sua derivada f’, o que pode ser uma tarefa não muito agradável, factor que deve ser tomado em conta na escolha deste método. Do lado positivo, temos que ele converge quadraticamente. 5.5.2. Erros e Convergência Nem sempre é convergente. Teorema 5.5.1. Se numa vizinhança (a,b) do zero z suficientemente pequena se verificar que f ε C2[a,b] e que para todo o ξ nessa vizinhança 0 < m1 <= |f’(ξ )| <= M1 e |f’’(ξ )| <= M2 então o método converge, e o erro satisfaz a relação: |ek+1| <= M|ek|2 com M = M2/2m1 Se o zero não for simples, a ordem de convergência do método de Newton degrada-se, sendo possível demonstrar que, no caso de zeros duplos, a convergência é apenas linear. Exemplo 5.5.1. (fazer) 5.6. MÉTODO DE MULLER Os métodos apresentados aproximavam a curva a uma recta. esta ideia pode ser generalizada no sentido de usar polinómios de grau superior para aproximar a função f. Os únicos candidatos são de segundo grau. O método de Muller consiste precisamente em interpolar a curva y=f(x) por uma parábola passando por 3 pontos e tomar a intersecção desta com o eixo dos xx como uma nova aproximação do zero. O polinómio interpolador é, como sabemos: p(x) = f(xk) + f[xk, xk-1](x-xk) + f[xk, xk-1, xk-2](x-xk)(x-xk-1) dado que (x – xk)(x – xk-1) = (x – xk)2 + (x-xk) xk – xk-1) a parábola interpoladora pode escrever-se na seguinte forma: p(x) = yk + ck (x – xk) + dk (x – xk) 2 em que, por simplicidade, se fez: yk = f(xk) dk = f[xk, xk-1, xk-2] ck = f[xk, xk-1] + dk(xk – xk-1) = f[xk, xk-1] + f[xk, xk-2] – f[xk-1, xk-2] Esta parabola intersecta o eixo dos xx nos pontos a determinados por: - ck + - (ck 2 – 4yk dk)1/2 a = xk + ----------------------------------------------------------------------2 dk uma forma equivalente, mas menos susceptível a erros de arredondamento, é: 2 yk a = xk - ------------------------------------------------------------------------ck + - (ck 2 – 4yk dk)1/2 Esta expressão fornece 2 valores, um correspondendo ao sinal +, e o outro, ao sinal -, pelo que se põe a questão de saber qual deles se deve tomar. É usual escolher o sinal que produza o maior  valor absoluto para o denominador no segundo membro desta expressão. Se nas iterações intermédias aparecerem valores complexos, despreza-se a parte imaginária de xk+1

A ordem de convergência é aprox. 1,84, quase idêntica à do método de Newton, com a vantegm adicional de não recorrer a derivadas. 5.7. UTILIZAÇÃO DE INTERPOLAÇÃO INVERSA x = g(y) Para aparecer alguma diferença, esta deve aparecer quando se interpola g por polinómios de segundo grau como no método de Muller. A parábola interpoladora de g nos pontos (yk-2, xk-2), (yk-1, xk-1) e (yk, xk) é: x = xk + g[yk, yk-1](y-yk) + g[yk, g[ yk, yk-1, yk-2](y – yk)(y – yk-1) A intersecção desta curva com o eixo dos xx obtém-se fazendo simplesmente y=0, donde se extrai que: xk+1 = xk – g[yk, yk-1]yk + g[yk, yk-1, yk-2]yk yk-1 Se tivermos em conta a expressão das diferenças divididas e desenvolvermos, chegamos a: yk ykyk-1 1 1 xk+1 = [ xk - ----------------------------------------------- ] + --------------- [ -------------- - ----------------- ] f[xk, xk-1] yk-2 – yk f[xk-1, xk-12] f[xk, xk-1] podemos ver que o primeiro parêntesis do segundo membro desta expressão corresponde exactamente ao método da secante. O segundo termo constituirá, portanto, uma correcção de segunda ordem, como, aliás, é de esperar. Nesta expressão, contrariamente ao que sucedia com o método de Muller, o valor de xk+1 é obtido sem qualquer ambiguidade e, se os valores iniciais forem reais, a sucessão {xk} é também real. p é também aprox. 1,84. 5.8. ITERAÇÃO DE PONTO FIXO 5.8.1. Descrição do Método Todos os métodos estudados até este momento se reportavam à equação não-linear f(x)=0. nesta secção vamos abordar a solução de equações não-lineares escritas desta feita na forma: x = g(x) (5.8.1.) Este modo não constitui qualquer restrição relativamente ao caso f(x)=0, pois é sempre possível transformar esta equação numa do tipo da 5.8.1., pondo, por ex, x = x + f(x) = g8x) A expressão 5.8.1. sugere imediatamente o seguinte esquema iterativo: xk+1 = g(xk), k=0,1,... Há casos em que o método é convergente e outros em que é divergente. O ponto z solução de (5.8.1.) é um ponto que a função g transforma nele próprio, ou, por outras palavras, um ponto que permanece fixo sob a transformação g. 5.8.2. Erros e Convergência Definição 5.8.1. Uma função g diz-se contractiva no intervalo I = [a,b] se existir uma constante M com 0<=M<1 tal que: |g(x1) – g(x2)| <= M |x1 – x2| para todos os x1, x2 ε I

CAP. 6 – SISTEMAS DE EQUAÇÕES LINEARES: LI NEARES: MÉTODOS DIRECTOS 6.1. Introdução x são as incógnitas, y são os coeficientes e b são os segundos membros do sistema de equações. Frequentemente é vantajoso colocar o sistema sob a forma matricial. Os métodos de solução de sistemas de equação lineares (sel) costumam ser classificados em 2 categorias: os métodos directos que permitem obter a solução de qualquer sistema com um nº finito de operações aritméticas e os métodos iterativos caso contrário. 6.1.1. Notação e Nomenclatura Definição 6.1.1. Uma matriz A de dimensão mxn é um quadro de números aij, com i01,...,m e j=1,...,n. Estes números são designados por elementos da matriz e costumam dispor-se em m linhas e n colunas. Uma matriz de dimensão x1 designa-se por vector coluna, e uma matriz de dimensão 1xn designa-se por vector linha. Se m=n a matriz diz-se quadrada; se m0!n diz-se rectangular. Os elementos da forma aii dizem-se elementos diagonais e forma a diagonal principal ou simplesmente diagonal da matriz A. Os elementos que se dispõem paralelamente à diagonal forma as codiagonais, que podem ser as supradiagonais (as que estão acima da diagonal) ou subdiagonais (abaixo). Definição 6.1.2. Uma matriz diz-se: Nula – se todos os elementos forem iguais a zero Trapezoidal Superior – se aij=0, qualquer i>j e inferior... Triangular Superior – se for simultaneamente quadrada e trapezoidal superior; triangular inferior  se... Diagonal se for simultaneamente triangular superior e inferior D=diag(d1,d2,...,dn). A matriz identidade I é uma matriz diagonal cujos elementos são todos iguais a 1 Uma matriz de ordem n diz-se simétrica se aij=aji i,j=1,2,...,n e anti-simétrica se aij=-aji i,j=1,2,...,n Positiva se aij>0 para todos os valores de i,j e este facto é denotado por A>0; negativa se... 6.1.2. Operações Com Matrizes Definição 6.1.3. A e B dizem-se iguais se tiverem as mesmas dimensões e os seus elementos homólogos forem iguais, isto é, aij=bij Definição 6.1.4. A soma é definida por cij=aij+bij (mesma dimensão para A e B ... e C) É comut comutati ativa va e asso associa ciativ tiva: a: A+B=B+ A+B=B+A A (A+B)+C=A+(B+C) Definição 6.1.5. O produto de um escalar α por uma matriz A é definida por bij=αaij. Goza das propriedades (αβ)A=α(βA) (α+β)A = αA+βA α(A+B) = αA + αB 1A = A 0A=0 As definições 4 e 5 tornam o conjunto das matrizes um espaço linear  Definição 6.1.6.  p a b Seja A ε Rm*p e B ε Rp*n. Então o produto de A por B é a matriz C=AB dada por cij = k =1 ik kj Para i=1,...,m, j=1,...,n Esta definição impõe que o nº de colunas de A seja igual ao nº de linhas de b. Matrizes com esta propriedade costumam designar-se por conformes à multiplicação. Prop oprrieda daddes da mult ultipli plicação (AB)C=A(BC) A(B+C) = AB + AC (A+B)C = AC + BC α(AB) = (αA)B Não é comutativa.

∑

Definição 6.1.7. A transposta de uma matriz A é uma matriz B cujos elementos são dados por: bij=aji e a notação usada é AT Goza das propriedades (AT)T = A (A+B)T = AT + BT (αA)T = αAT (AB)T BT AT Daqui resulta que uma matriz simétrica é igual à sua transposta e uma matriz anti-simétrica é igual à simétrica da sua transposta

Definição 6.1.8. Uma matriz A tal que ATA=AA T=I diz-se ortogonal Definição 6.1.9. (importante) Se existir uma matriz X tal que AX=XA = I , diz-se que X é a inversa da matriz A, a qual se denota por A-1 Por esta definição vemos que só as matrizes quadradas podem aspirar a possuir inversa. No próximo teorema compilaremos algumas propriedades da matriz inversa que usaremos muitas vezes.

Teorema 6.1.1. Seja A uma matriz de ordem n, e A-1 a sua inversa. Então são verdadeiras a seguintes proposições: A inversa, caso exista, é única (A-1)-1 = A (AT)-1 = (A-1)T Se a e B forem matrizes de ordem n invertíveis, então AB é invertível e (AB)-1 = B-1 A-1 Algumas classes de matrizes muito frequentes nas aplicações são introduzidas na próxima definição.

Definição 6.1.10. – Uma matriz A diz-se: definida positiva se xT Ax > 0 para qualquer X=!0 (x é um vector) definida semipositiva definida negativa e seminegativa Definição 6.1.11 Uma matriz A diz-se diagonal dominante por linhas se n

| aii | >=

∑| a

ij|

  j =1 =j ! i

para i de 1 até n

analogamente diz-se que A é diagonal dominante por colunas, bastando trocar aij por aji Se se tirar o igual dizem-se estritamente dominantes. Matrizes Particionadas As matrizes podem ser particionadas em submatrizes, ficando assim com os seus elementos a ser  também eles matrizes. As vantagens são a compactação da escrita e poder fazer-se operações sobre essas submatrizes. Um modo natural de particionar é tomar para submatrizes as respectivas linhas e colunas: Ai.T será a sua linha e a.j a sua coluna, o que torna a matriz representável por uma coluna e por  uma linha respectivamente. Quando se considera matriz A particionada por colunas, a posmultiplicação de A por um vector x adquire um novo e importante significado. De facto: x1 Ax = (a.1 a.2 ... a.n) x2 = a .1 x1 + a .2 x2 +…+ a .n xn … xn O que significa que Ax é a combinação linear das colunas de A cujos coeficientes são as componentes do vector x.

Identicamente a premultiplicação de x por A adquire um significado semelhante, ou seja, xTA é a combinação linear das linhas de A cujos coeficientes são as componentes do vector x. 6.1.3. Teoria dos Sistemas de Equações Lineares Um dos problemas centrais da Álgebra Linear consiste na pesquisa das condições que garantem a existência e unicidade das soluções de sistemas de equações do tipo Ax = b (6.1.2.) A formulação destas condições simplifica-se, e a sua interpretação adquire novos e importantes significados se identificarmos as linhas e as colunas de A como vectores e Ax como uma combinação linear das colunas de A em que as componentes xi do vector x desempenham o papel de coeficientes da combinação linear. Há 4 espaços lineares associados a uma matriz, que introduziremos. Definição 6.1.12. O espaço das colunas da matriz A, denotado por R(A), é o conjunto de vectores y ε Rm da forma y=Av, v ε Rn, isto é, R(A) é constituído por todas as combinações lineares das colunas de A. Analogamente para os espaço das linhas. Com esta definição o problema de resolver o sistema pode ser formulado de outra maneira: procurar a combinação linear das colunas de A que reproduz o segundo membro b. Daqui extraímos a conclusão que necessário e suficiente para a existência de soluções que b pertença a R(A). Definição 6.1.13. O espaço de nulidade ou núcleo de A, denotado por N(A), é o conjunto formado pelas soluções do sistema homogéneo Ax=0 O espaço de nulidade esquerdo resulta de que N(AT) = {v ε Rm : AT v = 0} Mas AT v = vT A = 0, com o vector v a aparecer do lado esquerdo de ª Os 4 conjuntos de vectores acabados de introduzir são espaços lineares, pelo que faz sentido falar da sua dimensão. Teorema 6.1.2. Seja r = r(A) o nº de colunas de A linearmente independentes. Então são verdadeiras as seguintes proposições: DimR(A) = r dimR(AT)=r dimN(A)=n-r dimN(AT)=m-r  É um resultado surpreendente que o espaço de colunas e de linhas acabem por ter a mesma dimensão, a que se dá o nome de característica da matriz A. Ao valor n-r costuma chamar-se nulidade de A ou deficiência da característica de A. Teorema 6.1.3. (Existência) O sistema 6.1.2. tem pelo menos uma solução x ε Rn para qualquer  segundo membro b ε Rm sse r=m<=n. (Unicidade) – a solução é única sse r=m=n Como se pode ver só as quadradas podem aspirar a ter solução única e daí irmos concentrar a nossa atenção nelas. Teorema 6.1.4. a solução do sistema 6.1.2. existe e é única sse o sistema homogéneo Ax=0 possuir solução trivial x=0 como única soluçãp, ou seja, sse a nulidade de A for zero. Teorema 6.15. O sistema 6.1.12 tem solução única sse qualquer das duas condições equivalentes for válida: A-1 existir  Det A =! 0 Teorema 6.1.6. (Regra de Cramer) – Seja A uma matriz invertível e Ai a matriz que se obtém de A substituindo nesta coluna i pelo segundo segundo membro b. Então, xi=detAi / detA com i de 1 a n. A simplicidade desta fórmula esconde o facto de que ela é totalmente inadequada como método numérico geral de solução. O cálculo de um determinante de ordem n requer (n-1)n! Multiplicações e n!-1 somas ou subtracções. 6.1.4. Sistemas Triangulares

Os métodos de solução que iremos estudar baseiam-se na ideia de transformar os sistema dado noutro, que seja mais fácil de resolver. Uma classe de sistemas fáceis de resolver são aqueles cujas matrizes são triangulares. Então temos de resolver Lx=b com L triangular inferior, de um modo geral as soluções podem ser  encontradas por substituições sucessivas (descendentes) ou ascendentes no caso de uma matriz triangular superior. Se o fizermos obtemos: i −1

Xi = (bi -

∑ l x ) / lii ij

j

com i de 1 a n

  j =1

A complexidade é de O(n2). Duas propriedades importantes das matrizes triangulares apresentam-se a seguir: Teorema 6.1.7. Seja T uma matriz triangular inferior (superior). Então, T é invertível sse os seus elementos diagonais forme diferentes de zero. Neste caso a inversa de T também é uma matriz triangular inferior (superior). 6.2. REDUÇÃO A SISTEMAS TRIANGULARES Os métodos directos tiram partido do facto de as matrizes triangulares serem fáceis de resolver, transformando o sistema original Ax=b num outro sistema A’x=b´, mas cuja matriz A´seja triangular. Definição 6.2.1. dois sistemas de equações lineares dizem-se equivalentes se possuírem o memso conjunto de soluções, ou seja, se os conjuntos {x} e [x} das soluções de Ax=b e A´x´=b ´forem iguais. Teorema 6.2.1. É condição suficiente para dois sistemas Ax e=b e A´x´=b´serem equivalentes que exista uma matriz C invertível tal que A´=CA e b´=Cb . se A for invertível, então a condição é também necessária. É fácil descobrir um conjunto de operações que não só garante a equivalência dos sistemas como também possibilita que a matriz transformada assuma a forma triangular pretendida. Definição 6.2.2. Dá-se a designação de operações elementares sobre as linhas de uma matriz às seguintes operações: Permutação de duas linhas Multiplicação de uma linha por um nº m diferente de zero Soma a uma linha do produto de outra linha por um nº m Podemos verificar através de exemplos que estas operações correspondem a premultiplicar a matriz A e o segundo membro b por certas matrizes invertíveis. Permutação de linhas Se premultiplicar a matriz identidade à qual se trocou a 1ª e 2ª linhas, por A, obtém-se A com as linhas 1 e 2 trocadas. Também a posmultiplicação produz os mesmos resultados. Multiplicação de uma linha por um nº m=!0 Se quiser multiplicar a linha 2 por m, basta premultiplicar (ou posmultiplicar) pela matriz

1  M=  0 0 

0

0  

  1    

m 0   0

Sendo M uma matriz diagonal, a condição M=!0 garante a sua invertibilidade. Soma a uma linha do produto de outra linha por um nº m Ex. Para somar à linha 2 de uma matriz A de ordem 3 o produto da linha 1 por m, podemos premultiplicar A por 

1 m  0 

0 0  

  1    

1 0   0

Matrizes desta forma vão ser muito úteis adiante, pelo que vamos proceder a uma generalização. Definição 6.2.4. Uma matriz triangular inferior elementar de ordem n e índice k é uma matriz da forma M = I – meTk em que m é um vector cujas componentes satisfazem mi=0, mi=0, i=1,...,k e ek é n o k-ésimo vector da base canónica de R Uma outra designação designação atribuída atribuída a estas matrizes matrizes é a de transformação transformação de Gauss. O seu tipo é o seguinte (ex. Para n=6 e k=3)

1 0  0  0 0  0

0

0

0 0

0  

1

0

0 0

0  

0

1

0 0

0  

0

−m4 −m5 −m6

1

0  

0 0

0

0 1

   

0  

 

0 0 1    

As (estas) matrizes elementares triangulares inferiores são invertíveis e o cálculo da sua inversa é muito simples, como iremos ver. Teorema 6.2.2. Seja M uma matriz triangular inferior elementar. Então M é invertível e a sua inversa é M-1 = I + meTk Podemos concluir que todas as operações elementares sobre linhas da matriz correspondem a premultiplicar A por matrizes invertíveis e, portanto, conduzem a sistemas equivalentes. Em segu seguid idaa vamo vamoss mo most stra rarr qu quee um enc encad adeam eamen ento to ade adequ quado ado des desta tass ope opera raçõ ções es pe perm rmititee transformar A numa triangular superior. 6.3. MÉTODO DE GAUSS Designaremos por A(k) a matriz transformada de A na etapa k. A primeira operação a efectuar consiste em anular o elemento a21, o que, como vimos, se pode fazer multiplicando multiplicando a primeira linha por m21 = a21/a11 e subtraindo o resultado à segunda linha. Este processo designa-se por condensação (do elemento a21). A segunda linha de A e de b modificam-se de acordo com as seguintes expressões: a2j(2) = a2j – m21a1j j=2,...,n b2(2) = b2-m21b1 Este processo pode ser repetido para todos os elementos da coluna 1 abaixo da diagonal, utilizando para a condensação do elemento ai1 o multiplicador mi1 = ai1 / a11 O resultado de todas estas condensações produz a matriz transformada A(2) (cuja primeira coluna abaixo da diagonal é formada por zeros) e produz ainda b(2) Podemos agora proceder com a matriz obtida (2) do mesmo modo que com (1) condensando a segunda coluna abaixo da diagonal e assim sucessivamente. As expressões gerais obtidas são: mik=aik(k) / akk(k) | (k+1) (k) (k) aij = aij – mik akj , j=k+1,…,n | I = k+1,…,n (6.3.1.) (k+1) (k) (k) bi = bi – mik bk | Os elementos akk são conhecidos por elementos pivot e a respectiva linha por linha pivot. A submatriz de A(k) constituída pelas suas n-k últimas linhas e n-k últimas colunas é designada por  submatriz activa. Deste modo, a condensação de Gauss consiste na subtracção de um múltiplo apropriado da linha pivot às linhas que lhe estão abaixo. Depois de ter A(n) = A´ é só utilizar o método das substituições ascendentes para chegar Às soluções. Do ponto de vista da programação há a realçar que os elementos aij da matriz original deixam de ser necessários assim que sofrem a sua primeira transformação, pelo que as matrizes A (k) podem ir sendo gravadas em cima da original para poupar memória.

Algoritmo 6.3.1. (Condensação de Gauss) para k01,...,n-1 fazer: para i=k+1,...,n fazer: mik = aik/akk para j = k+1,...,n fazer: aij = aij – mik*akj fim do ciclo j bi = bi – mik + bk fim do ciclo i fim do ciclo k A complexidade é de O(2n3/3) o que é muito menor do que a regra de Cramer. Se os elementos pivot não forem zero, mediante uma alteração simples, podemos viabilizar o método de Gauss. Vários Segundos Membros 6.4. FACTORIZAÇÕES TRIANGULARES 6.4.1. Factorização LU A condensação de gauss consiste numa sequência de operações elementares sobre as linhas das matrizes A(k), as quais, como vimo, podem ser expressas sob a forma de premultiplicações desta matriz por matrizes triangulares elementares. Por ex. A(2) = M(1) A(1) em que:

 1  −m 21 M(1) =   M   −m1

0 L

0  

1 L

0  

M O

0 L

 

    1   M

= I – m(1) e1T

Como podemos concluir do método de Gauss, em geral: A(k+1) = M(k) A(k) M(k) = I – m(k)ekT com k=1,…,n-1 (k) em que m é o vector constituído pelos multiplicadores e cujas componentes mi(k), i=1,...,k são nulas. No final teremos que: _  (n) A = A = M (n-1) M(n-2)...M(1) A Mas o produto de matrizes triangulares inferiores de diagonal unitária é também uma matriz triangular inferior de diagonal unitária. Pelo que MA = A barra. Do teorema 6.1.1. sabemos que M(k) é invertível, logo M também o é e dá: M-1 = (M(n-1)...M(1))-1 = (M(1))-1...(M(n-1))-1 = (I+m(1)e1T)...(I+m(n-1)en-1T) = I + m(1)e1T + ... + m(n-1)en-1T O 2º membro é também uma matriz triangular inferior de diagonal unitária cujos elementos no triângulo inferior são os multiplicadores mij surgidos na condensação de Gauss. Pondo, por ser  mais sugestivo:  _  A = U e M-1 = L ficamos com A = LU em que L é triangular inferior e U triangular superior – é a factorização triangular.

Exemplo 6.4.1. Factorizar a seguinte matriz pelo método de Gauss

A=

1 2  0 

0 2 1

2  

  1    

1   = A(1)

calculando os multiplicadores e os elementos transformados pelas expressões 6.3.1.:

mik(k)=aik(k) / akk(k) | (k+1) (k) (k) (k) aij = aij – mik akj , j=k+ j=k+1, 1,…, …,n n | i = k+1, k+1,…, …,n n (6.3 (6.3.1 .1.) .) (k+1) (k) (k) bi = bi – mik bk |

e colocando cada multiplicador na mesma posição do elemento da matriz A cuja condensação provocou, temos sucessivamente as seguintes matrizes transformadas: - multiplicadores:

m21(1) = a21(1) / a11(1) = 2/1 = 2 m31 = a31(1) / a11(1) = 0/1 = 0 a22(2) = a22(1) – m21(1) * a12(1) = 2 – (2*0) = 2 a23(2) = a23(1) – m21(1) * a13(1) = 1 – (2*2) = -3 a32(2) = a32(1) – m31(1) * a12(1) = 1 – (0*0) = 1 a33(2) = a33(1) – m31(1) * a13(1) = 1 – (0*2) = 1

A=

1 2  &&&  0

0 2 1

2  

(2) −3    = A 1    

Com os multiplicadores sublinhados.

m32(3) = a32(2) / a22(2) = ½ a33(3) = a33(2) – m32(2) * a23(2) = 1 – (1/2 * -3) 5/2 1  A = &2&& 0 

0 2 1/ 2

  −3     = A(3) 5 / 2    2

1  Então achámos L = A-1 =  2 0   1 0 2      e U =  0 2 −3    0 0 5 / 2     

0 1 1/ 2

0  

  1     0  

Do ponto de vista da resolução do sistema de equações: Ax=b (LU)x=b L(Ux) = b introduzindo o vector auxiliar y temos que

Ly=b e Ux=y

Conforme se vê, a determinação da solução x, uma vez obtida a factorização LU, passa pela solução em sucessão de dois sistemas triangulares: um inferior e outro superior.

Esta observação permite aliás revelar uma vantagem do método de factorização sobre o método de Gauss. A factorização envolve apenas a matriz A e não o segundo membro b, intervindo este exclusivamente na fase de solução dos sistemas triangulares. isto significa que, uma vez factorizada A, podemos resolver tantos sistemas quantos quisermos, à custa apenas de substituições ascendentes e descendentes. Como a factorização é a aprte mais cara em termos de nº de operações, há vantagem relativamente ao métodode Gauss aplicado integralmente a cada segundo membro.

... 6.4.4. FACTORIZAÇÃO LDU Uma pequena modificação das factorizações de Doolittle e de Crout, que tem a vantagem de atribuir às matrizes triangulares L e U um papel totalmente idêntico, é considerar que esta são ambas matrizes triangulares com diagonal unitária e tomar a factorização da seguinte forma: A=LDU em que D é uma matriz diagonal. O método de Doolittle consiste em absorver a matriz D e em U, ou seja, considerr a factorização A=L(DU) enquanto o de Crout é A=(LD)U 6.6. Cálculo da Inversa e do Determinante 6.6.1. Cálculo da Matriz Inversa Resolvendo, por qualquer dos métodos, o seguinte sistema: AX=I A inversa da matriz só deve ser calculada se explicitamente pedido e não como forma de resolver  sistemas de equações lineares.

6.6.2. Cálculo do Determinante Uma vez obtida a factorização de A, o cálculo do seu determinante pode fazer-se com um custo adicional desprezável. Suponhamos que havíamos obtido a factorização: PA=LDU Então, recordando que o determinante do produto é igual ao produto dos determinantes: (detP)(detA) = (detL)(detD)(detU) Mas, tendo em atenção a propriedade 6.2.2. das matrizes de permutação, e que o determinante de matrizes de diagonal unitária é igual a 1, podemos escrever que: n

det A = (det P) (det D) = (det P) ∏ d ii i =1

O valor de det P é +1 se o nº de permutações de linhas levadas a cabo tiver sido par e –1, se tiver  sido ímpar. Mesmo com matrizes de dimensão moderada pode dar overflow. Então deve-se incluir  no programa uma forma automática de escalar os valores de forma a que não ultrapassem esses limites. Como sabemos, uma matriz é singular sse det A=0. No entanto é incorrecto testar a singularidade com base no valor do determinante. O problema com o determinante é de que ele é afectado pela escala dos elementos da matriz enquanto a singularidade da matriz não o é. Para testar a singularidade então deve-se (ver 6.5.1.) tentar encontrar um pivot diferente de zero. Se nãofor possível, é singular. Mas os erros de arredondamento podem atrapalhar. Um indicador  que é correntemente utilizado é o do valor absoluto do pivot |akk(k)|. Um valor baixo, da ordem da unidade de arredondamento, é sinal de que a matriz é singular ouestá muito próximo de o ser, pelo que há que proceder com a devida cautela, nomeadamente no que diz respeito à fiabilidade dos resultados.