Exemplos selecionados de cancelamento catastrófico

Use o botão abaixo para reportar erros ou dar sugestões.

Cálculo Numérico - Versão GNU Octave

2.9 Exemplos selecionados de cancelamento catastrófico

Exemplo 2.9.1. Considere o seguinte processo iterativo: $\begin{array}{rcl} x^{(1)} & = & \frac{1}{3} & (2.156) \\ x^{(n + 1)} & = & 4 x^{(n)} - 1, n = 1, 2, \dots & (2.157) \end{array}$

Observe que $x^{(1)} = \frac{1}{3}$ , $x^{(2)} = 4 \cdot \frac{1}{3} - 1 = \frac{1}{3}$ , $x^{(3)} = \frac{1}{3}$ , ou seja, temos uma sequência constante igual a $\frac{1}{3}$ . No entanto, ao calcularmos no computador, usando o sistema de numeração double, a sequência obtida não é constante e, de fato, diverge. Faça o teste no GNU Octave, colocando:

>> x = 1/3

e itere algumas vezes a linha de comando:

>> x = 4*x-1

Para compreender o que acontece, devemos levar em consideração que o número $\frac{1}{3} = 0, \bar{3}$ possui um representação infinita tanto na base decimal quanto na base binária. Logo, sua representação de máquina inclui um erro de arredondamento. Seja $ϵ$ a diferença entre o valor exato de $\frac{1}{3}$ e sua representação de máquina, isto é, ${\tilde{x}}^{(1)} = \frac{1}{3} + ϵ$ . A sequência efetivamente calculada no computador é: $\begin{array}{rcl} {\tilde{x}}^{(1)} & = & \frac{1}{3} + ϵ & (2.158) \\ {\tilde{x}}^{(2)} & = & 4 x^{(1)} - 1 = 4 (\frac{1}{3} + ϵ) - 1 = \frac{1}{3} + 4 ϵ & (2.159) \\ {\tilde{x}}^{(3)} & = & 4 x^{(2)} - 1 = 4 (\frac{1}{3} + 4 ϵ) - 1 = \frac{1}{3} + 4^{2} ϵ & (2.160) \\ ⋮ & (2.161) \\ {\tilde{x}}^{(n)} & = & \frac{1}{3} + 4^{(n - 1)} ϵ & (2.162) \end{array}$

Portanto o limite da sequência diverge,

lim_{x \to \infty} | {\tilde{x}}^{(n)} | = \infty

(2.163)

Qual o número de condicionamento desse problema?

Exemplo 2.9.2. Observe a seguinte identidade

f (x) = \frac{(1 + x) - 1}{x} = 1

(2.164)

Calcule o valor da expressão à esquerda para $x = 1 0^{- 12}$ , $x = 1 0^{- 13}$ , $x = 1 0^{- 14}$ , $x = 1 0^{- 15}$ , $x = 1 0^{- 16}$ e $x = 1 0^{- 17}$ . Observe que quando $x$ se aproxima do $ϵ$ de máquina a expressão perde o significado. Veja a Figura 2.1 com o gráfico de $f (x)$ em escala logarítmica.

Figura 2.1: Gráfico na função do Exemplo 2.9.2.

Figura 2.2: Gráfico de

{(1 + \frac{1}{n})}^{n}

em função de

n

em escala linear-logarítmica variando de

1 0^{0}

até

1 0^{18}

. Veja o Exemplo 2.9.3.

Exemplo 2.9.3. Neste exemplo, estamos interessados em compreender mais detalhadamente o comportamento da expressão

{(1 + \frac{1}{n})}^{n}

(2.165)

quando $n$ é um número grande ao computá-la em sistemas de numeral de ponto flutuante com acurácia finita. Um resultado bem conhecido do cálculo nos diz que o limite de (2.165) quando $n$ tende a infinito é o número de Euler:

lim_{n \to \infty} {(1 + \frac{1}{n})}^{n} = e = 2, 718281828459 . . .

(2.166)

Sabemos também que a sequência produzida por (2.165) é crescente, isto é:

{(1 + \frac{1}{1})}^{1} < {(1 + \frac{1}{2})}^{2} < {(1 + \frac{1}{3})}^{3} < \dots

(2.167)

No entanto, quando calculamos essa expressão no Octave, nos defrontamos com o seguinte resultado:







$n$	${(1 + \frac{1}{n})}^{n}$	$n$	${(1 + \frac{1}{n})}^{n}$











$1$	$2, 0000000000000$	$1 0^{2}$	$2, 7048138294215$
$2$	$2, 2500000000000$	$1 0^{4}$	$2, 7181459268249$
$3$	$2, 3703703703704$	$1 0^{6}$	$2, 7182804690957$
$4$	$2, 4414062500000$	$1 0^{8}$	$2, 7182817983391$
$5$	$2, 4883200000000$	$1 0^{10}$	$2, 7182820532348$
$6$	$2, 5216263717421$	$1 0^{12}$	$2, 7185234960372$
$7$	$2, 5464996970407$	$1 0^{14}$	$2, 7161100340870$
$8$	$2, 5657845139503$	$1 0^{16}$	$1, 0000000000000$
$9$	$2, 5811747917132$	$1 0^{18}$	$1, 0000000000000$
$10$	$2, 5937424601000$	$1 0^{20}$	$1, 0000000000000$

Podemos resumir esses dados no gráfico de ${(1 + \frac{1}{n})}^{n}$ em função de $n$ , veja a Figura 2.9.

Observe que quando $n$ se torna grande, da ordem de $1 0^{15}$ , o gráfico da função deixa de ser crescente e apresenta oscilações. Observe também que a expressão se torna identicamente igual a $1$ depois de um certo limiar. Tais fenômenos não são intrínsecos da função $f (n) = {(1 + 1 ∕ n)}^{n}$ , mas oriundas de erros de arredondamento, isto é, são resultados numéricos espúrios. A fim de pôr o comportamento numérico de tal expressão, apresentamos abaixo o gráfico da mesma função, porém restrito à região entre $1 0^{14}$ e $1 0^{16}$ .

Para compreendermos melhor por que existe um limiar $N$ que, quando atingido torna a expressão do exemplo acima identicamente igual a $1$ , observamos a sequência de operações realizadas pelo computador:

n \to 1 ∕ n \to 1 + 1 ∕ n \to {(1 + 1 ∕ n)}^{n}

(2.168)

Devido ao limite de precisão da representação de números em ponto flutuante, existe um menor número representável que é maior do que 1. Este número é $1 +$ eps, onde eps é chamado de épsilon de máquina e é o menor número que somado a 1 produz um resultado superior a 1 no sistema de numeração usado. O épsilon de máquina no sistema de numeração double vale aproximadamente $2, 22 \times 1 0^{- 16}$ . No GNU Octave, o épsilon de máquina é a constante eps. Observe que:

>> printf(’%1.25f\n’, 1+eps)
1.0000000000000002220446049

Quando somamos a $1$ um número positivo inferior ao épsilon de máquina, obtemos o número 1. Dessa forma, o resultado obtido pela operação de ponto flutuante $1 + n$ para $0 < n < 2, 22 \times 1 0^{- 16}$ é 1.

Portanto, quando realizamos a sequência de operações dada em (2.168), toda informação contida no número $n$ é perdida na soma com $1$ quando $1 ∕ n$ é menor que o épsilon de máquina, o que ocorre quando $n > 5 \times 1 0^{15}$ . Assim, $(1 + 1 ∕ n)$ é aproximado para $1$ e a última operação se resume a $1^{n}$ , o que é igual a $1$ mesmo quando $n$ é grande.

Um erro comum é acreditar que o perda de significância se deve ao fato de $1 ∕ n$ ser muito pequeno para ser representado e é aproximando para $0$ . Isto é falso, o sistema de ponto de flutuante permite representar números de magnitude muito inferior ao épsilon de máquina. O problema surge da limitação no tamanho da mantissa. Observe como a seguinte sequência de operações não perde significância para números positivos x muito menores que o épsilon de máquina:

n \to 1 ∕ n \to 1 ∕ (1 ∕ n)

(2.169)

compare o desempenho numérico desta sequência de operações para valores pequenos de $n$ com o da seguinte sequência:

n \to 1 + n \to (1 + n) - 1 .

(2.170)

Finalmente, notamos que quando tentamos calcular ${(1 + \frac{1}{n})}^{n}$ para $n$ grande, existe perda de significância no cálculo de $1 + 1 ∕ n$ . Para entendermos isso melhor, vejamos o que acontece no GNU Octave quando $n = 7 \times 1 0^{13}$ :

>> format(’long’)
>> n=7e13
n =  70000000000000
>> 1/n
ans =    1.42857142857143e-14
>> y=1+1/n
y =  1.00000000000001

Observe a perda de informação ao deslocar a mantissa de $1 ∕ n$ . Para evidenciar o fenômenos, observamos o que acontece quando tentamos recalcular $n$ subtraindo $1$ de $1 + 1 ∕ n$ e invertendo o resultado:

>> y-1
ans = 1.42108547152020e-14
>> 1/(y-1)
ans = 70368744177664

Exemplo 2.9.4 (Analogia da balança). Observe a seguinte comparação interessante que pode ser feita para ilustrar os sistemas de numeração com ponto fixo e flutuante: o sistema de ponto fixo é como uma balança cujas marcas estão igualmente espaçadas; o sistema de ponto flutuante é como uma balança cuja distância entre as marcas é proporcional à massa medida. Assim, podemos ter uma balança de ponto fixo cujas marcas estão sempre distanciadas de $100$ g ( $100$ g, $200$ g, $300$ g, ..., $1$ Kg, $1, 1$ Kg,...) e outra balança de ponto flutuante cujas marcas estão distanciadas sempre de aproximadamente um décimo do valor lido ( $100$ g, $110$ g, $121$ g, $133$ g, ..., $1$ Kg, $1, 1$ Kg, $1, 21$ Kg, ...) A balança de ponto fixo apresenta uma resolução baixa para pequenas medidas, porém uma resolução alta para grandes medidas. A balança de ponto flutuante distribui a resolução de forma proporcional ao longo da escala.

Seguindo nesta analogia, o fenômeno de perda de significância pode ser interpretado como a seguir: imagine que você deseje obter o peso de um gato (aproximadamente $4$ Kg). Dois processos estão disponíveis: colocar o gato diretamente na balança ou medir seu peso com o gato e, depois, sem o gato. Na balança de ponto flutuante, a incerteza associada à medida do peso do gato (sozinho) é aproximadamente $10 %$ de $4$ Kg, isto é, $400$ g. Já a incerteza associada à medida da uma pessoa (aproximadamente $70$ Kg) com o gato é de $10 %$ do peso total, isto é, aproximadamente $7$ Kg. Esta incerteza é da mesma ordem de grandeza da medida a ser realizada, tornado o processo impossível de ser realizado, já que teríamos uma incerteza da ordem de $14$ Kg (devido à dupla medição) sobre uma grandeza de $4$ Kg.

Exercícios resolvidos

ER 2.9.1. Deseja-se medir a concentração de dois diferentes oxidantes no ar. Três sensores eletroquímicos estão disponíveis para a medida e apresentam a seguintes respostas:

v_{1} = 270 [A] + 30 [B], v_{2} = 140 [A] + 20 [B] e v_{3} = 15 [A] + 200 [B]

(2.171)

as tensões $v_{1}$ , $v_{2}$ e $v_{3}$ são dadas em $m V$ e as concentrações em $m i l i m o l ∕ l$ .

a)

Encontre uma expressão para os valores de

[A]

[B]

em termos de

v_{1}

v_{2}

e, depois, em termos de

v_{1}

v_{3}

. Dica: Se

a d \neq b c

, então a matriz

A

dada por

A = [\begin{matrix} a & b \\ c & d \end{matrix}]

(2.172)

é inversível e sua inversa é dada por

A^{- 1} = \frac{1}{a d - b c} [\begin{matrix} d & - b \\ - c & a \end{matrix}] .

(2.173)

b)

Sabendo que incerteza relativa associada às sensibilidades dos sensores 1 e 2 é de

2 %

e que a incerteza relativa associada às sensibilidades do sensor 3 é

10 %

, verifique a incerteza associada à medida feita com o par

1 - 2

e o par

1 - 3

. Use

[A] = [B] = 10 m i l i m o l ∕ l

. Dica: Você deve diferenciar as grandezas

[A]

[B]

em relação aos valores das tensões.

Solução. Em ambos casos, temos a seguinte estrutura:

[\begin{matrix} S_{11} & S_{12} \\ S_{21} & S_{22} \end{matrix}] [\begin{matrix} [A] \\ [B] \end{matrix}] = [\begin{matrix} v_{1} \\ v_{2} \end{matrix}]

(2.174)

De forma que

[\begin{matrix} [A] \\ [B] \end{matrix}] = {[\begin{matrix} S_{11} & S_{12} \\ S_{21} & S_{22} \end{matrix}]}^{- 1} [\begin{matrix} v_{1} \\ v_{2} \end{matrix}] = \frac{1}{S_{11} S_{22} - S_{12} S_{21}} [\begin{matrix} S_{22} & - S_{12} \\ - S_{21} & S_{11} \end{matrix}] [\begin{matrix} v_{1} \\ v_{2} \end{matrix}]

(2.175)

Portanto $\begin{array}{rcl} [A] & = & \frac{S_{22} v_{1} - S_{12} v_{2}}{S_{11} S_{22} - S_{12} S_{21}} & (2.176) \\ [B] & = & \frac{- S_{21} v_{1} + S_{11} v_{2}}{S_{11} S_{22} - S_{12} S_{21}} & (2.177) \end{array}$

Usando derivação logarítmica, temos $\begin{array}{rcl} \frac{1}{[A]} \frac{\partial [A]}{\partial S_{11}} & = & - \frac{S_{22}}{S_{11} S_{22} - S_{12} S_{21}} & (2.178) \\ \frac{1}{[A]} \frac{\partial [A]}{\partial S_{12}} & = & - \frac{v_{2}}{S_{22} v_{1} - S_{12} v_{2}} + \frac{S_{21}}{S_{11} S_{22} - S_{12} S_{21}} = - \frac{[A]}{[B]} \cdot \frac{S_{22}}{S_{11} S_{22} - S_{12} S_{21}} & (2.179) \\ \frac{1}{[A]} \frac{\partial [A]}{\partial S_{21}} & = & \frac{S_{12}}{S_{11} S_{22} - S_{12} S_{21}} & (2.180) \\ \frac{1}{[A]} \frac{\partial [A]}{\partial S_{22}} & = & \frac{v_{1}}{S_{22} v_{1} - S_{12} v_{2}} - \frac{S_{11}}{S_{11} S_{22} - S_{12} S_{21}} = \frac{[A]}{[B]} \cdot \frac{S_{12}}{S_{11} S_{22} - S_{12} S_{21}} & (2.181) \end{array}$

e $\begin{array}{rcl} \frac{1}{[B]} \frac{\partial [B]}{\partial S_{11}} & = & \frac{v_{2}}{- S_{21} v_{1} + S_{11} v_{2}} - \frac{S_{22}}{S_{11} S_{22} - S_{12} S_{21}} = \frac{[B]}{[A]} \frac{S_{21}}{S_{11} S_{22} - S_{12} S_{21}} & (2.182) \\ \frac{1}{[B]} \frac{\partial [B]}{\partial S_{12}} & = & \frac{S_{21}}{S_{11} S_{22} - S_{12} S_{21}} & (2.183) \\ \frac{1}{[B]} \frac{\partial [B]}{\partial S_{21}} & = & - \frac{v_{1}}{- S_{21} v_{1} + S_{11} v_{2}} + \frac{S_{21}}{S_{11} S_{22} - S_{12} S_{21}} = - \frac{[B]}{[A]} \frac{S_{11}}{S_{11} S_{22} - S_{12} S_{21}} & (2.184) \\ \frac{1}{[B]} \frac{\partial [B]}{\partial S_{22}} & = & - \frac{S_{11}}{S_{11} S_{22} - S_{12} S_{21}} & (2.185) \\ (2.186) \end{array}$

E o erro associado às medidas pode ser aproximado por $\begin{array}{rcl} \frac{1}{[A]} δ_{[A]} & = & |\frac{1}{[A]} \frac{\partial [A]}{\partial S_{11}}| δ_{S_{11}} + |\frac{1}{[A]} \frac{\partial [A]}{\partial S_{12}}| δ_{S_{12}} + |\frac{1}{[A]} \frac{\partial [A]}{\partial S_{21}}| δ_{S_{21}} + |\frac{1}{[A]} \frac{\partial [A]}{\partial S_{22}}| δ_{S_{22}} & (2.187) \\ = & \frac{1}{|det S|} [S_{22} δ_{S_{11}} + \frac{[A]}{[B]} S_{22} δ_{S_{12}} + S_{12} δ_{S_{21}} + \frac{[A]}{[B]} S_{12} δ_{S_{22}}] & (2.188) \end{array}$

Analogamente, temos: $\begin{array}{rcl} \frac{1}{[B]} δ_{[B]} & = & \frac{1}{|det S|} [\frac{[B]}{[A]} S_{21} δ_{S_{11}} + S_{21} δ_{S_{11}} + \frac{[B]}{[A]} S_{11} δ_{S_{21}} + S_{11} δ_{S_{22}}] & (2.189) \end{array}$

onde não se indicou $| S_{i j} |$ nem $| [.] |$ pois são todos positivos.

Fazemos agora a aplicação numérica:
Caso do par 1-2:

det S = |\begin{matrix} 270 & 30 \\ 140 & 20 \end{matrix}| = 1200

(2.190)

\begin{array}{rcl} \frac{1}{[A]} δ_{[A]} & = & \frac{1}{1200} [20 \times 270 \times 2 % + 20 \times 30 \times 2 % + 30 \times 140 \times 2 % + 30 \times 20 \times 2 %] & (2.191) \\ = & \frac{216}{1200} = 0.18 = 18 % & (2.192) \\ \frac{1}{[B]} δ_{[B]} & = & \frac{1}{1200} [140 \times 270 \times 2 % + 140 \times 30 \times 2 % + 270 \times 140 \times 2 % + 270 \times 20 \times 2 %] & (2.193) \\ = & \frac{426}{1200} = 0.355 = 35.5 % & (2.194) \end{array}

Caso do par 1-3:

det S = |\begin{matrix} 270 & 30 \\ 15 & 200 \end{matrix}| = 53550

(2.195)

\begin{array}{rcl} \frac{1}{[A]} δ_{[A]} & = & \frac{1}{53550} [200 \times 270 \times 2 % + 200 \times 30 \times 2 % + 30 \times 15 \times 10 % + 30 \times 200 \times 10 %] & (2.196) \\ = & \frac{1804, 6}{52550} \approx 0.0337 = 3.37 % & (2.197) \\ \frac{1}{[B]} δ_{[B]} & = & \frac{1}{53550} [15 \times 270 \times 2 % + 15 \times 30 \times 2 % + 270 \times 15 \times 10 % + 270 \times 200 \times 10 %] & (2.198) \\ = & \frac{5895}{53550} \approx 0.11 = 11 % & (2.199) \end{array}

Conclusão, apesar de o sensor $3$ apresentar uma incerteza cinco vezes maior na sensibilidade, a escolha do sensor $3$ para fazer par ao sensor $1$ parece mais adequada.

$♢$

Exercícios

E 2.9.1. Considere as expressões:

\frac{exp (1 ∕ μ)}{1 + exp (1 ∕ μ)}

(2.200)

\frac{1}{exp (- 1 ∕ μ) + 1}

(2.201)

com $μ > 0$ . Verifique que elas são idênticas como funções reais. Teste no computador cada uma delas para $μ = 0, 1$ , $μ = 0, 01$ e $μ = 0, 001$ . Qual dessas expressões é mais adequada quando $μ$ é um número pequeno? Por quê?

Resposta. Quando

μ

é pequeno,

e^{1 ∕ μ}

é um número grande. A primeira expressão produz um ”overflow” (número maior que o máximo representável) quando

μ

é pequeno. A segunda expressão, no entanto, reproduz o limite

1

quando

μ \to 0 +

E 2.9.2. Encontre expressões alternativas para calcular o valor das seguintes funções quando $x$ é próximo de zero.

a): $f (x) = \frac{1 - cos (x)}{x^{2}}$
b): $g (x) = \sqrt{1 + x} - 1$
c): $h (x) = \sqrt{x + 1 0^{6}} - 1 0^{3}$
d): $i (x) = \sqrt{1 + e^{x}} - \sqrt{2}$ Dica: Faça $y = e^{x} - 1$

Resposta. a)

\frac{1}{2} + \frac{x^{2}}{4!} + O (x^{4})

; b)

x ∕ 2 + O (x^{2})

; c)

5 \cdot 1 0^{- 4} x + O (x^{2})

; d)

\frac{\sqrt{2}}{4} y + O (y^{2}) = \frac{\sqrt{2}}{4} x + O (x^{2})

E 2.9.3. Use uma identidade trigonométrica adequada para mostrar que:

\frac{1 - cos (x)}{x^{2}} = \frac{1}{2} {(\frac{sen (x ∕ 2)}{x ∕ 2})}^{2} .

(2.202)

Analise o desempenho destas duas expressões no computador quando $x$ vale $1 0^{- 5}$ , $1 0^{- 6}$ , $1 0^{- 7}$ , $1 0^{- 8}$ , $1 0^{- 9}$ , $1 0^{- 200}$ e $0$ . Discuta o resultado. Dica: Para $| x | < 1 0^{- 5}$ , $f (x)$ pode ser aproximada por $1 ∕ 2 - x^{2} ∕ 24$ com erro de truncamento inferior a $1 0^{- 22}$ .

Resposta. A expressão da direita se comporta melhor devido à retirada do cancelamento catastrófico em

x

em torno de

0

E 2.9.4. Reescreva as expressões:

\sqrt{e^{2 x} + 1} - e^{x} e \sqrt{e^{2 x} + x^{2}} - e^{x}

(2.203)

de modo que seja possível calcular seus valores para $x = 100$ utilizando a aritmética de ponto flutuante ("Double") no computador.

Resposta. Possíveis soluções são:

\begin{array}{rcl} \sqrt{e^{2 x} + 1} - e^{x} & = & \sqrt{e^{2 x} + 1} - e^{x} \cdot \frac{\sqrt{e^{2 x} + 1} + e^{x}}{\sqrt{e^{2 x} + 1} + e^{x}} & (2.204) \\ = & \frac{e^{2 x} + 1 - e^{2 x}}{\sqrt{e^{2 x} + 1} + e^{x}} = \frac{1}{\sqrt{e^{2 x} + 1} + e^{x}} & (2.205) \end{array}

e, de forma análoga: $\begin{array}{rcl} \sqrt{e^{2 x} + x^{2}} - e^{x} = \frac{x^{2}}{\sqrt{e^{2 x} + x^{2}} + e^{x}} . & (2.206) \end{array}$

E 2.9.5. Na teoria da relatividade restrita, a energia cinética de uma partícula e sua velocidade se relacionam pela seguinte fórmula:

E = m c^{2} (\frac{1}{\sqrt{1 - {(v ∕ c)}^{2}}} - 1),

(2.207)

onde $E$ é a energia cinética da partícula, $m$ é a massa de repouso, $v$ o módulo da velocidade e $c$ a velocidade da luz no vácuo dada por $c = 299792458 m ∕ s$ . Considere que a massa de repouso $m = 9, 10938291 \times 1 0^{- 31} K g$ do elétron seja conhecida com erro relativo de $1 0^{- 9}$ . Qual é o valor da energia e o erro relativo associado a essa grandeza quando $v = 0, 1 c$ , $v = 0, 5 c$ , $v = 0, 99 c$ e $v = 0, 999 c$ sendo que a incerteza relativa na medida da velocidade é $1 0^{- 5}$ ?