Cálculo Num érico – UNICAMPbiloti/an/precisaofinita.pdf · do sistema numérico; t, quantos...

Introdução Cancelamento Equação quadrática Somas infinitas Norma Euclidiana

Computação em precisão finita

Ricardo [email protected]

Cálculo Numérico – UNICAMP

1S/2021

http://bit.ly/blt-tj2oD

http://bit.ly/blt-tj2oD Ricardo Biloti Computação em precisão finita


Licença

Este trabalho é licenciado sob os termos da Licença InternacionalCreative Commons Atribuição-NãoComercial-CompartilhaIgual 4.0.

Para ver uma cópia desta licença, visitehttp://creativecommons.org/licenses/by-nc-sa/4.0/.


Seus direitos e deveres são:

• Você é livre para copiar e redistribuir este material, em qualquer meio ou formato,para adaptá-lo, transformá-lo ou utilizá-lo para construir seu próprio material.

• Você deve dar os créditos apropriados, fornecendo link para a licença e indicando sealterações foram feitas. Você pode fazer isto de qualquer forma razoável, porém semtentar passar a ideia ou sugerir que o autor endosse suas alterações ou seu uso domaterial.

• Você não pode utilizar este material para fins comerciais.

• Se você alterar, transformar ou construir seu próprio material com base nestetrabalho, você deverá distribúı-lo sob a mesma licença usada no original.


Erros v https://youtu.be/ygKZLToNg9w

Se x é uma quantidade numérica e x̂ sua aproximação, então

Eabs(x̂) = |x − x̂ |

Erel (x̂) =|x − x̂ ||x | , (x 6= 0)

Eadm(x̂) =|x − x̂ |

L

onde L é uma dimensão caracteŕıstica do problema.


Repare que no cálculo do erro absoluto e do erro relativo sempre é necessário conhecer x .Ou seja, para medir se uma aproximação é boa ou não é necessário compará-la com o valorexato. Porém, em geral procuramos aproximações quando não conhecemos o valor exato.Como então medir o erro de uma aproximação? Em problemas práticos isto nem sempre éde fato posśıvel.

Vários métodos numéricos contudo fornecem estimativas para o erro absoluto e/ou relativodas aproximações que produzem.

Quando se deparar com um método numérico ou com uma aproximação numérica, você devesempre se perguntar se é posśıvel fornecer ou se é conhecida uma estimativa para o erro. Docontrário, como julgar a qualidade da aproximação?

Também é usual (e mais viável) adimensionalizar o erro absoluto utilizando para isso umadimensão caracteŕıstica do problema. Por exemplo, se o problema for estimar a distânciaentre duas cidades próximas, uma dimensão caracteŕıstica razoável seria L = 1 km. Se oproblema for estimar a população de um páıs, L = 1 milhão de habitantes faz sentido. Jápara trabalhar com a altura de pessoas L poderia ser escolhido como 1.70 m.

https://youtu.be/ygKZLToNg9w


Exemplo

I x = 1.00000, e x̂ = 1.00499

Eabs(x̂) = 4.99 · 10−3 Erel (x̂) = 4.99 · 10−3

I x = 9.00000, e x̂ = 8.99501

Eabs(x̂) = 4.99 · 10−3 Erel (x̂) = 5.54 · 10−4


Um erro absoluto de mesma ordem pode ser mais ou menos significativo dependendo dasgrandezas envolvidas. Por exemplo, para quem quer comprar uma toalha de mesa, um metroa menos ou a mais faz muita diferença, mas para quem está interessado na distância entreduas cidades, um erro na casa dos metros é despreźıvel.

No exemplo dado, apesar dos erros absolutos serem da mesma ordem (≈ 10−3), os errosrelativos são de ordens diferentes. Isto permite concluir que no segudo caso a aproximação émelhor que no primeiro caso.


Tipos de erros v https://youtu.be/ygKZLToNg9w

I Erros de medição/aquisição

I Erros de representação

I Erros de cálculo em precisão finita

I Erros introduzidos por algoritmos numéricos


• Erros de medição são inerentes à aquisição de dados experimentais, podendos serapenas minorados, mas nunca evitados. São causados por falha humana, deequipamento, imprecisões do experimento, etc.

• Erros de representação são igualmente inevitáveis. Eles surgem sempre que os dadosdo problema são digitalizados. Dentre todas as fontes de erros, essa é a menosproblemática. Em geral, erros de representação são muitas ordens de magnitudemenor que os erros originados por outras causas.

• Erros introduzidos pelo cálculo em precisão finita são o objeto de estudo deste tópicodo curso. Estudaremos qual o impacto na qualidade das quantidades calculadascausado pelo fato das contas serem feitas em no computador.

• Por fim, há erros introduzidos pelo emprego de métodos numéricos que apenasaproximam a solução de um determinado problema. Esta fonte de erro será analizadano decorrer do curso, sempre que um novo método for abordado.

https://youtu.be/ygKZLToNg9w


Representação em ponto flutuante

x =√

3 = 1.732050807568877...

Num sistema de ponto flutuante (SPF)

x̂ = fl(x) = 1.732051= 0.1732051 · 10+01

I Base: 10I Mantissa: 0.1732051I Expoente: 01


As quatro caracteŕısticas que definem um sistema de ponto flutuante (ou SPF) são: β, a basedo sistema numérico; t, quantos d́ıgitos significativos, ou mantissa, são armazenados; L, qualo menor expoente representável, e U, qual o maior expoente representável. Simbolicamente,representamos um sistema de ponto flutuante por F(β, t, L,U).

Para evitar múltiplas representações para o mesmo número, convenciona-se que expoenteseja definido de tal forma que o primeiro d́ıgito de um número de ponto flutuante sejasempre zero e o segundo d́ıgito seja sempre diferente de zero. Do contrário o número 1.3por exemplo podeŕıa tanto ser representado como 0.13000 · 101 como 0.01300 · 102 ou como0.00130 · 103, e assim por diante.

Por exemplo, num sistema de ponto flutuante, de base 10, capaz de armazenar 7 d́ıgitos paraa mantissa e expoentes entre −99 e 99, fl(

√3) = 1.732051 seria representado como

fl(√

2) = 0.1732051 · 10+01,

onde fl(x) é a notação para a representação de ponto flutuante do número x .

No computadores atuais, o usual é que o sistema de ponto flutuante utilize base 2, armazene52 bits (52 d́ıgitos binários), e aceite expoentes entre −1022 e 1023.


Exemplos

Num sistema de base 10, 5 d́ıgitos para mantissa, e expoente entre−99 e 99

π = 3.14159265358979... fl(π) = 3.1416c = 299 792.458 km/s fl(c) = 299 790.000 km/s

Qual o menor e o maior número representável em módulo?

0.10000 · 10−99 e 0.99999 · 1099


Nos exemplos do slide, como π e c seriam representados no sistema de ponto flutuante, ouseja, qual seria a mantissa e o expoente para cada um?

No sistema de ponto flutuante descrito, qual seria a menor e a maior distância entre doisnúmeros consecutivos representáveis?

No Octave, o menor e o maior múmeros representáveis (positivos) podem ser conferidos comos comandos realmin e realmax.


Precisão × acurácia

PrecisãoErro cometido em operações algébricas elementares

AcuráciaErro presente em quantidade aproximadas

https://youtu.be/hRAFPdDppzs?list=PL2jCIH8XciGksBisuY7I5PnVtv1NV56VR


É importante distinguir precisão de acurácia. Enquanto que o primeiro termo se refere a umapropriedade do sistema de ponto flutuante da máquina, o segundo termo diz respeito à todaestratégia utilizada para a cálculo da quantidade aproximada.

Dizemos que uma máquina é muito precisa se os erros em operações como somas/subtraçõese produtos/divisões for pequeno. Já o adjetivo acurado se presta a qualificar o resultado finalobtido pelo algoritmo. Um algoritmo para o cálculo de tensões em uma estrutura metálicapor exemplo envolve milhões de operações algébricas elementares. Sua qualidade não édeterminada apenas pela qualidade com que estas operações são executadas.

Por fim, cabe destacar que estes dois termos são empregados de maneira confusa e sãomuitas vezes intercambiados. O comum é que a palavra precisão seja utilizada nos doiscontextos.

https://youtu.be/hRAFPdDppzs?list=PL2jCIH8XciGksBisuY7I5PnVtv1NV56VRhttps://youtu.be/hRAFPdDppzs?list=PL2jCIH8XciGksBisuY7I5PnVtv1NV56VR


Precisão

A precisão de um SPF é quantificada pela unidade dearredondamento u. Esse é o menor número u tal que

fl(x � y) = (x � y)(1 + δ), |δ| < u,

onde � representa +, −, × ou ÷.

Isso é o mesmo que dizer que∣∣∣∣fl(x � y)− (x � y)(x � y)∣∣∣∣ = |δ| < u,

u é o maior erro relativo que podemos esperar em uma operaçãoaritmética.


A precisão de um sistema de ponto flutuante é quantificada pela unidade de arredondamentou. Esse número limita qual o erro relativo máximo que pode ocorrer em uma operaçãoalgébrica elementar.

O padrão IEEE-754-1985 (revisado em 2008), Standard for Binary Floating-Point Arithmetic,é aplamente adotado por fabricantes de processadores. Este padrão normatiza o sistema deponto flutuante implementado em processadores.

A diferença entre precisão simples e precisão dupla é o espaço utilizado para o armazenamentodo número de ponto flutuante. Em precisão simples são utilizados 32 bits, enquanto que 64bits são necessários para a representação em precisão dupla. O esquema abaixo descrevecomo esses bits são distribúıdos para representar a mantissa (f ) e o expoente (e), além dosinal (s) de um número de ponto flutuante.

1 8 23Precisão simples: s e f

1 11 52Precisão dupla: s e f

A unidade de arredondamento, assim como outras quantidades como os maiores e menoresnúmeros representáveis, e os resultados esperados de ações como arredondamentos, trunca-mentos, etc, são propriedades do sistema de ponto flutuante.

Steve Hollash escreveu um bom texto introdutório de sobre o padrão IEEE-754-1985(http://steve.hollasch.net/cgindex/coding/ieeefloat.html).


Estimando a unidade de arredondamento

A unidade de arredondamento é o menor número positivo de pontoflutuante u tal que

fl(1 + u) > 1

No padrão IEEEI u ≈ 1.19 · 10−07 (em precisão simples)I u ≈ 2.22 · 10−16 (em precisão dupla)


Uma forma de determinar a unidade de arrendondamento da máquina é descobrir qual omenor valor para u, representável no SPF, tal que fl(1 + u) ainda é maior que 1.


Precisão 6⇒ acurácia

Alta precisão não é suficiente para garantir resultados acurados


O importante é ter em mente que máquinas muito precisas não bastam para obter resultadosacurados.

Pense, por exemplo, que não basta um ter bisturi a laser (instrumento de alta precisão)para que uma cirurgia seja bem sucedida (resultado acurado). Por outro lado, sem um bombisturi, dificilmente uma cirurgia seria bem sucedida.

Via de regra, as máquinas atuais são bem precisas. Mesmo assim, para muitos problemas édif́ıcil conseguir boas soluções aproximadas.


Erro de cancelamento v https://youtu.be/xSoGZYWYMaw

Sistema de ponto flutuante com cinco d́ıgitos significativos.

Queremos calcular

49213 + 31.728− 49244 = 0.728


Para efeito de exposição, vamos considerar um sistema de ponto flutuante simples (SPFS)capaz de representar apenas os cinco d́ıgitos mais significativos de um número e um únicod́ıgito para expoente. Neste SPFS, a unidade de arredondamento é u = 5 · 10−5. Paraver isso, perceba que se fl(1+5·10−5) = fl(1.00005) = 1.0001, por conta do arredondamento.

Vamos analisar a seguinte conta simples:

49213 + 31.728− 49244 = 0.728

Repare que todos os número envolvidos são representados corretamente no nosso sistema deponto flutuante simplificado.

https://youtu.be/xSoGZYWYMaw


Erro de cancelamento

49213 + 31.728− 49244 = 0.728

49213+ 31.728

49244.728

fl(49213 + 31.728) = 49245, Erel = 6 · 10−6

fl(49245− 49244) = 1, Erel = 0


As operações são realizadas da esquerda para a direita. A primeira soma é feita e o resultadoé um número com oito d́ıgitos significativos. Sendo assim, apenas os cinco mais significativossão preservados em nosso sistema de ponto flutuante.

A conta é conclúıda com a subtração de 49244 de 49245, cujo resultado exato 1 é corretamenteobtido.


Erro de cancelamento

x = 49213 + 31.728− 49244 = 0.728

x̂ = fl(49213 + 31.728− 49244) = 1

Erel (x̂) = 0.374� 5 · 10−5 = u


Neste exemplo o valor correto x deveria ser 0.728, enquanto que o valor obtido foi 1. O errorelativo nesta aproximação é 0.374, bem maior que o erro embutido nas operações algébricaselementares executadas.

A única conta realizada com erro foi a adição inicial. Porém o erro relativo naquela operaçãofoi pequeno

|49245− 49244.728|49244.728

= 5.52 · 10−6,

o que é compat́ıvel com o sistema de ponto flutuante simples. Entretanto o erro relativo daaproximação como um todo é bem maior:

|1− 0.728|0.728

= 3.74 · 10−1.

Como isto é posśıvel? Onde este erro grosseiro foi originado?


Perda de d́ıgitos

Três d́ıgitos significativos perdidos

49213+ 31.728

49244.728

fl(49213 + 31.728) = 49213 + 32 = 49245

A subtração foi exata

fl(49245− 49244) = 1


O problema surgiu na adição inicial. O resultado numérico daquela operação só pôdepreservar os cinco d́ıgitos mais significativos dos oito d́ıgitos que compunham o resultadoexato. Isso não é um problema para a adição em si. De fato, o erro relativo nesta operaçãofoi de 5.52 · 10−6, menor que a unidade de arredondamento do SPFS. Entretanto aquelestrês d́ıgitos perdidos, que não eram importantes para o resultado da adição, passaram a serimportantes quando a subtração seguinte foi executada.

Na subtração, todos aqueles cinco d́ıgitos mais significativos foram cancelados. Os trêsd́ıgitos perdidos seriam então novamente necessários, porém não há mais como recuperá-los.

Neste sentido, a subtração apenas evidenciou um problema que foi gerado pela perda ded́ıgitos significativos em um passo anterior.


Fórmula de Heron para a área de triângulos

A =√

s(s − a)(s − b)(s − c), s = (a + b + c)/2

10-15

10-13

10-11

10-9

1 10-1 10-2 10-3 10-4 10-5 10-6 10-7Erro

rela

tivo

no c

álcu

lo d

a ár

ea

Altura do triângulo


Heron, no ano 60, escreveu a área de um triângulo em termos do comprimentos de seuslados. Numericamente, a área computada desta forma torna-se mais imprecisa a medida queo triângulo torna-se mais achatado. Para tais triângulos, s será aproximadamente igual a umdos lados, digamos a. Nesse caso, a fórmula tal como apresentada, sofre pelo cancelamentode d́ıgitos significativos na subtração de s − a.

No gráfico, foi exibido o erro relativo no cálculo da área de um triângulo retângulo comhipotenusa 4 e altura progressivamente menor.

Uma maneira de evitar esse problema é computar a área pela fórmula

A = 14√

(a + (b + c))(c − (a − b))(c + (a − b))(a + (b − c))

onde os lados do triângulo foram ordenados pelo comprimento, isto é, a ≥ b ≥ c. Para evitaro cancelamento de d́ıgitos, o cálculo deve ser feito na ordem indicada pelos parênteses.


Exerćıcio

Para que valores de x , as expressões abaixo podem sofrer por errosde cancelamento?

Reescreva as expressões abaixo de maneira a reduzir posśıveis erros.

1.√

x2 + 1− x2.√

1 + x − 13. (1− cos x)/ sin x4.√

(1− cos x)/2


A idéia para este exerćıcio é procurar uma forma de rescrever as expressões de maneira aevitar o cancelamento de d́ıgitos.

Por exemplo, no caso de√

x2 + 1 − x , se x for muito grande teremos que√

x2 + 1 ≈ x ,e portanto a subtração destas quantidades implicará no cancelamento de d́ıgitos significativos.

Esta expressão pode porém ser reescrita como√x2 + 1− x =

[√x2 + 1− x

] √x2 + 1 + x√

x2 + 1 + x= 1√

x2 + 1 + x,

que não sofre de cancelamento uma vez que não há a subtração de quantidades próximas.

Por exemplo, se x = 109, fl(√

x2 + 1− x)

= 0, enquanto que fl(

1/[√

x2 + 1 + x ])

=5 · 10−10.


Equação quadrática

ProblemaEncontrar as duas ráızes reais de

x2 − bx + c = 0

Se b2 − 4c ≥ 0,

r = b ±√

b2 − 4c2

ExerćıcioEscreva um algoŕıtmo para computar as ráızes reais de umaequação quadrática.


Considere o problema de encontrar as ráızes de uma equação quadrática. Como existefórmula fechada para as ráızes desta equação, basta utilizá-la. Apenas por conveniência,considere a equação normalizada de maneira que o termo quadrático tenha coeficiente 1.

Seu algoritmo não precisa ser escrito em nenhuma linguagem espećıfica. Basta que fique claroquais operações devem ser realizadas e em que ordem.


Equação quadrática

Vamos aplicar seu algoritmo para este exemplo.

x2 − bx + c = 0

b = 4.7379100021 c = 0.0016199351

r1 = 4.7375680682... r2 = 0.0003419340...


Para esse exemplo, com os coeficientes b e c dados, r1 e r2 são as duas ráızes desta equação(todos os d́ıgitos exibidos estão corretos).


Sequência de cálculo

r = b −√

b2 − 4c2

b = 4.7379100021 c = 0.0016199351

1. fl(b2) = 2.2448 · 10+12. fl(4c) = 6.4797 · 10−33. fl(b2 − 4c) = 2.2442 · 10+14. fl(

√b2 − 4c) = 4.7373 · 10+0

5. fl(b −√

b2 − 4c) = 6.0000 · 10−46. fl((b −

√b2 − 4c)/2) = 3.0000 · 10−4


Vamos acompanhar a sequência de cálculo necessária para o compto da menor ráız. Todasas operações são realizadas no nosso sistema de ponto flutuante simplificado.


Erro

Erel (r̂2) =|3.0000 · 10−4 − 3.4193 · 10−4|

3.4193 · 10−4 = 1.2263·10−1 � 10−4 = u


O erro relativo é bem superior à unidade de arredondamento do SPFS (que é uma medidapara o erro relativo máximo esperado em operações elementares).



r = b −√

b2 − 4c2

b = 4.7379100021 c = 0.0016199351

1. fl(b2) = 2.2448 · 10+12. fl(4c) = 6.4797 · 10−33. fl(b2 − 4c) = 2.2442 · 10+14. fl(

√b2 − 4c) = 4.7373 · 10+0

5. fl(b −√

b2 − 4c) = 6.0000 · 10−46. fl((b −

√b2 − 4c)/2) = 3.0000 · 10−4




r = b −√

b2 − 4c2

b = 4.7379100021 c = 0.0015

1. fl(b2) = 2.2448 · 10+12. fl(4c) = 6.0000 · 10−33. fl(b2 − 4c) = 2.2442 · 10+14. fl(

√b2 − 4c) = 4.7373 · 10+0

5. fl(b −√

b2 − 4c) = 6.0000 · 10−46. fl((b −

√b2 − 4c)/2) = 3.0000 · 10−4


De fato, repare que o valor encontrado para a ráız seria o mesmo se, ao invés de utilizar1.6199 ·10−3 para c, tivéssemos utilizado c = 1.5 ·10−3. Ou seja, quatro d́ıgitos significativosde c foram perdidos no decorrer das operações. Esta perda só foi sentida no passo 5, quandoa subtração cancelou os d́ıgitos mais significativos e aqueles anteriormente perdidos passariama ser novamente importantes.


Maior raiz

r = b +√

b2 − 4c2

b = 4.7379100021 c = 0.0016199351

1. fl(b2) = 2.2448 · 10+12. fl(4c) = 6.4797 · 10−33. fl(b2 − 4c) = 2.2442 · 10+14. fl(

√b2 − 4c) = 4.7373 · 10+0

5. fl(b +√

b2 − 4c) = 9.4752 · 10+06. fl((b +

√b2 − 4c)/2) = 4.7376 · 10+0


Se ao invés de calcular a menor ráız, calculássemos a maior, o passo 5 seria uma adição aoinvés de uma subtração, e não haveria mais a necessidade de manter os d́ıgitos perdidos nopasso 3.


Erro

Erel (r̂1) =|4.7376− 4.7375680682|

4.7375680682 = 6.74 · 10−6


Veja que o erro relativo no compto da maior ráız é bem menor, e compat́ıvel com o erro dearredondamento do SPFS.

Como fazer para estimar então a menor ráız?


Segunda raiz

r2 =cr1

= 3.4193 · 10−4 Erel (r̂2) = 1.1600 · 10−5


Um alternativa inteligente é trocar o algoritmo para o cálculo desta ráız. Utilizando a relaçãor1r2 = c, podemos calcular a menor ráız sem fazer qualquer subtração que evidenciaria umaperda de d́ıgitos significativos anterior.


Estratégia

r1 =b + sign(b)

√b2 − 4c

2 r2 =cr1

ProblemasI b2 ≈ 4cI overflow ou underflow em b2


A melhor estratégia então seria sempre calcular primeiro a ráız que não tem problema comcancelamento de d́ıgitos. A outra ráız seria calculada através da relação r1r2 = c.

Isto é o melhor que pode ser feito com esta expressão para as ráızes de uma equaçãoquadrática. Entretanto isto ainda não resolve todos os problemas.

Se b2 ≈ 4c ainda haverá cancelamento de d́ıgitos significativos que não pode ser evitado.Por fim, o cálculo de b2 ainda pode apresentar overflow ou underflow.


Exemplo: Somas infinitas v https://youtu.be/ShC0-DiyXiY

Considere a somaS =

∞∑k=1

1k2 =

π2

6

Computacionalmente

S ≈ SN =N∑

k=1

1k2

|S − SN | diminui a medida que N aumenta.


Computacionalmente, não podemos somar infinitos termos. Por isto, fixamos um certo N erealizamos a soma até este ı́ndice.

Como a série é convergente e como seus termos são todos positivos, sabemos que |S − SN |vai a zero, monotonicamente.

Logo, quão maior for N melhor é a aproximação de S por SN . Como cient́ısta numéricos,podemos fazer um experimento para observar isto.

https://youtu.be/ShC0-DiyXiY


Algoritmo

I s ← 1

I k ← 2

I enquanto k ≤ N,I s ← s + 1/k2

I k ← k + 1

https://youtu.be/6hfOvs8pY1k?list=PL2jCIH8XciGksBisuY7I5PnVtv1NV56VR


Uma forma de validar ou refutar a tese proposta é conduzir um experimento numérico bemcontrolado.

O algoritmo mais simples e intuitivo para realizar computacionalmente esta soma acumulaparcela por parcela os termos da soma.

O śımbolo ← significa uma atribuição, ou seja, k ← k + 1, significa que à variável k seráatribúıdo o valor que ela tem atualmente acrescido de 1.

https://youtu.be/6hfOvs8pY1k?list=PL2jCIH8XciGksBisuY7I5PnVtv1NV56VRhttps://youtu.be/6hfOvs8pY1k?list=PL2jCIH8XciGksBisuY7I5PnVtv1NV56VR


Resultado

1

1e-09

1e-07

1e-05

1e-03

1e-01

0 5 10 15 20 25 30

|SN

- (π2

/6)|

log2(N)


Neste gráfico exibimos o erro∣∣∣SN − π26 ∣∣∣, para N = 2n, calculado pelo algoritmo proposto.

Podemos ver claramente que a partir de um certo valor de N, temos uma estagnação naredução do erro.

Para este produzir este gráfico as contas foram feitas em precisão simples. Isso nãocompromete a análise que estamos fazendo. Como veremos a seguir, apenas antecipa umfenônemo que também aconteceria se tivéssemos usado precisão dupla.

Perceba ainda que o menor erro obtido foi da ordem de 10−4, o que é bem superior àunidade de arredondamento em precisão simples (u ≈ 10−7).

A saber, para forçar que uma variável do Octave seja de precisão simples, inicialize-a comosingle. Por exemplo, a = single(0). Use o comando whos, para verificar o tipo de cadavariável.


Análise

s = S4095 = 1.64472532 . . .

Para k = 4096,

S4096 = fl(s + 1/40962) = fl(s + 2−24)

= fl(

s(

1 + 2−24

s

))

= s fl(

1 + 2−24

s

)= s

2−24s ≈ 3.5804 · 10

−8 < 1.1920 · 10−7 ≈ u


Depois de somadas 4095 parcelas, temos que s = S4095 = 1.64472532 . . .. A soma parcialseguinte, S4096, é calculada como

S4096 = fl(S4095 + 1/40962).

Como s já é um número de ponto flutuante, fl(s) = s e como fl(s ·p) = fl(s) fl(p), temos que

S4096 = s fl(

1 + 2−24

s

).

Porém, como 2−24/s é menor que a unidade de arredondamento (em precisão simples), oresultado de ponto flutuante da soma desta quantidade com 1 é 1, e portanto a soma nãose altera: S4096 = S4095.

O problema aqui é que o termo 1/40962 é muito pequeno em comparação com S4095 eportanto a sua soma é despreźıvel dentro do sistema de ponto flutuante.


Perda de d́ıgitos

A perda de d́ıgitos ocorre quando somamos números de grandezasmuito distintas.

SoluçãoI k ← N, s ← 0

I enquanto k ≥ 1,I s ← s + 1/k2

I k ← k − 1


A alternativa é acumular as parcelas na ordem inversa. Desta forma, cada nova parcela (maior)sempre será somada a um valor acumulado também crescente, não havendo a disparidade degrandezas que acontece no algoritmo ingênuo.


Comparação

1

1e-09

1e-07

1e-05

1e-03

1e-01

0 5 10 15 20 25 30

|SN

- (π2

/6)|

log2(N)

Ordem crescenteOrdem decrescente


Neste gráfico exibimos o erro∣∣∣SN − π26 ∣∣∣, para N = 2n, calculado pelo algoritmo crescente

(ingênuo) e calculado pelo algoritmo decrescente (ordem inversa). Enquanto que o erro sereduz apenas até N = 212 = 4096 no algoritmo ingênuo, no algoritmo decrescente é posśıvelchegar perto da precisão da máquina. Note que as contas foram feitas em precisão simples.


Conclusão

Perda de d́ıgitos significativos pode ocorrer sem que hajacancelamento.


Verificamos que a tese inicial, de que sempre que ao se acrescentar mais termos a soma, oresultado computacional torna-se mais preciso não era correta. De fato, a pergunta estavamal posta, no sentido de que computacionalmente, não se pode falar sobre uma fórmula, masapenas sobre a implementação da fórmula — um algoritmo.


Norma de Euclidiana

Se x = (x1, x2, . . . , xn)T ,

‖x‖2 =√

x21 + x22 + · · ·+ x2nAlgoritmoI k ← 1, s ← 0

I enquanto k ≤ n,I s ← s + (xk · xk )

I k ← k + 1

I s ←√

s


Outro exemplo de operação simples mas que deve ser feita com cuidado é o cálculo da normade um vetor.

A fórmula da norma euclidiana é facilmente traduzida em um algoritmo simples.


Exemplo 1

Em um sistema de ponto flutuante com cinco d́ıgitos significativose expoente entre −9 e 9,

x = (6 · 104, 8 · 104), ‖x‖2 = 105,

s = 3.6 · 109 + (8 · 104 × 8 · 104) (overflow)


Nesse exemplo simples, o cálculo da norma pelo algoritmo ingênuo apresentado, para um ve-tor com apenas duas componentes, já resulta em overflow. Repare que tantos as coordenadasdo vetor como o valor da norma são quantidades representadas no sistema de ponto flutuante.

Veja que o problema acontece quando 8 · 104 deve ser elevado ao quadrado, pois esse termosozinho já gera o overflow.


Exemplo 2

Em um sistema de ponto flutuante com cinco d́ıgitos significativose expoente entre −9 e 9,

x = (7 · 104, 6 · 104, 5 · 104), ‖x‖2 = 1.0488 · 105,

s = (7 · 104 × 7 · 104) + (6 · 104 × 6 · 104) + (5 · 104 × 5 · 104)(overflow)


Nesse exemplo, novamente todas as componentes do vetor podem ser representadas no sis-tema de ponto flutuante, seu valores ao quadrado também são representáveis, assim como ovalor da norma do vetor. Porém o acúmulo dos valores ao quadrado das componentes gera ooverflow.


Cálculo da norma

ProblemasI Overflow/underflow quando s = s + x2k

I Overflow/underflow no cálculo de x2k


Com os dois exemplos, vimos que tanto o passo de elevar uma componente ao quadrado,quanto o passo de acumular esses valores podem gerar overflow (assim como underflow).


Escalamento

x = (6 · 104, 8 · 104)

‖x‖2 =√

(6 · 104)2 + (8 · 104)2

=

√√√√(8 · 104)2 [(6 · 1048 · 104)2

+ 1]

= 8 · 104√

(6/8)2 + 1

= 105


A solução, no caso desse algoritmo, é escalar o vetor pela sua compomente de maior valorabsoluto antes de computar sua norma.


Algoritmo

I k ← 1, s ← 0, γ ← max |xk |,

I enquanto k ≤ n,I s ← s + (|xk |/γ) · (|xk |/γ)

I k ← k + 1

I s ← γ√

s


Essa versão do algoritmo não sofre dos problemas dos exemplos anteriores, sendo semprecapaz de computar a norma de um vetor, desde que o valor da norma seja representável nosistemas de ponto flutuante.

Porém, há um incoveniente. Esse algoritmo precisa percorrer o vetor duas vezes. A primeiradelas, apenas para poder descobrir o maior valor absoluto das componentes do vetor, de modoa realizar o escalamento.


Algoritmo

I k ← 1, γ ← 1, s ← 0

I enquanto k ≤ n,I Se γ > |xk |, então

I s ← s + (|xk |/γ) · (|xk |/γ)senão dnrm2 (BLAS-1)I s ← 1 + s(γ/|xk |) · (γ/|xk |)

I γ ← |xk |

I k ← k + 1

I s ← γ√

s


Essa última versão faz o escalamento sem precisar pecorrer o vetor duas vezes, mas simdescobrindo o fator de escala durante o processo e corrigindo-o, se necessário.

Esse é o algoritmo que de fato é utilizado para computar norma de vetores em bibliote-cas de algoritmos numéricos de qualidade, como a BLAS (do inglês, Basic Linear AlgebraSubprograms).

IntroduçãoCancelamentoEquação quadráticaSomas infinitasNorma Euclidiana

Cálculo Num érico – UNICAMPbiloti/an/precisaofinita.pdf · do sistema numérico; t, quantos...

Documents

Transcript of Cálculo Num érico – UNICAMPbiloti/an/precisaofinita.pdf · do sistema numérico; t, quantos...