Item a)
Inicialmente, é interessante apresentar os conceitos de precisão simples (ou única) e precisão dupla (ou single precision e double precision). Para isso, primeiro devemos responder a questão: como computadores armazenam números reais? Sabemos que toda informação é armazenada de forma binária: zeros ou uns. Então, a princípio, uma simples conversão de base poderia parecer adequada. Mas como representar números irracionais, por exemplo? Ou frações exatas na base decimal, mas que geram sequências infinitas na base binária? (como será visto mais a frente). Como padronizar a quantidade de memória que será alocada para armazenar cada número? (além de outras questões, como otimizar as operações aritméticas realizadas pelo processador.)
Como solução, adota-se o chamado sistema de ponto flutuante (floating-point system). Nele, cada número real \(N\) é representado por uma mantissa (\(m\)), um expoente (\(e\)), uma base (\(b\)), e (\(s\)), que define o sinal no número real, tal que:
\[
N=(-1)^s \times m \times b^e
\]
Em sistemas digitais, \(b=2\), e a mantissa e o expoente são obviamente armazenados na forma binária. Por exemplo, para representarmos o número real 5.75:
Conversão da base decimal para a base binária: \(5.75_{10} = 1\times2^2+0\times2^1+1\times2^0+1\times2^{-1}+1\times2^{-2} = 101.11_{2}\)
"Normalizando": \(101.11_{2} = 1.0111_{2}\times2^2\)
Ou seja, uma possível representação de \(5.75_{10}\) em um sistema de ponto flutuante binário é:
\[
5.75_{10} = (-1)^0\times1.0111_{2}\times2^{{10}_2}
\]
Em que \(s = 0_{2}, m = 1.0111_{2}, e = 10_{2} (=2_{10})\).
E em casos de números com uma quantidade bem maior de dígitos, ou mesmo infinita, como definir quantos bits devemos alocar para cada um desses elementos?
Aqui entra o conceito de precisão simples e precisão dupla. O Institute of Electrical and Electronics Engineers, por meio da norma IEEE 754, definiu requisitos mínimos para esses dois formatos, sendo que o formato de precisão simples utiliza 32 bits para armazenar um número real, enquanto o formato de precisão dupla, 64 bits. E a quantidade de bits para cada elemento da representação de ponto flutuante é divida da seguinte forma em cada padrão:
\[
\begin{array}{r|ccc|c}
Precisão & Sinal (s) & Expoente (e) & Mantissa (m) & Total(bits)\\
\hline
Simples & 1 & 8 & 23 & 32\\
Dupla & 1 & 11 & 52 & 64
\end{array}
\]
Existem outros detalhes da implementação desses formatos, como rotação do expoente para permitir expoentes negativos, e omissão do primeiro bit da mantissa. Esses detalhes podem ser conferidos nas referências ao final dessa resposta. Também é importante salientar que a IEEE 754 define requisitos mínimos, logo algumas implementações podem, internamente, utilizar precisões maiores, mas armazenar os resultados finais na precisão indicada.
Como curiosidade, no caso do formato de precisão simples, você pode testar como funciona a conversão de números reais para a representação de ponto flutuante de forma rápida em:
www.h-schmidt.net/FloatConverter/IEEE754.html
Note, por exemplo, que não é possível armazenar de forma precisa um valor real como \(0.1_{10}\) (base decimal) em uma representação de base binária, incluindo a de ponto flutuante. No caso de precisão simples, o valor real armazenado será:
0.100000001490116119384765625
ou seja, com erro na ordem de \(10^{-9}\). Isso ocorre pois \(0.1_{10}\) é representado por uma sequência infinita na forma \(0.00011001100110011...\) na base binária. Como temos uma limitação de bits para representar o número, haverá esse erro devido ao truncamento do valor até o limite de bits disponível. O problema é o mesmo da fração \(1/3=0.333...\) na base decimal: precisamos de uma sequência infinita de \(3\)'s para representá-la.
E \(0.1_{10}\), em precisão dupla, na realidade será armazenado como:
0.1000000000000000055511151231257827021181583404541015625
ou seja, um erro na ordem de \(10^{-18}\). Para verificar esse valor, no Python, basta executar:
from decimal import Decimal
Decimal.from_float(0.1)
Mas então por que, no Python, ao realizarmos uma operação como \(1/10\) no console, o resultado padrão exibido é exatamente 0.1? Pois a exibição do resultado geralmente é arredondada até certas casas decimais. Mas, para verificar, por exemplo, o problema que apontamos acima com relação as representações binárias, tente executar \(.1 + .1 + .1 == .3\) em um console Python, ou verifique o resultado da operação \(0.8 * 12\) no console padrão do Spyder, e surpreenda-se com a resposta...
Apesar do erro ser pequeno nessas representações, no caso de um grande número de operações complexas ou consecutivas, estes erros podem ser propagados e amplificados, gerando, em alguns casos, resultados finais com imprecisões em um nível indesejado. E assim, sem entrar em mais detalhes sobre as representações e operações aritméticas envolvendo ponto flutuante, já pode-se notar que a escolha de um ou de outro formato, simples ou duplo, pode determinar a precisão do resultado de operações envolvendo números reais.
Para o problema em questão, se resolvido por qualquer técnica não computacional em que carregamos as frações, quando existentes, até o final, encontramos como solução o vetor \({x} = [-1, 2, -3, 4]\).
Porém, para explorar as questões envolvendo precisão de ponto flutuante, utilizamos a implementação em Python descrita a seguir, que é um solver simples de sistemas lineares, em que dados \({A}\), \({b}\) e uma determinada precisão, retorna a solução \({x}\).
É importante salientar que, para simular precisão simples, alimentamos todas as operações algébricas do numpy com dados em precisão simples ("float32"). Dessa forma, o resultado de cada operação também será em precisão simples. Na codificação, em cada operação, ainda utilizamos o parâmetro dtype=precisao, mas não é necessário para nosso objetivo, visto que o numpy realizará a operação na precisão desejada caso os operandos estejam na mesma precisão.
Observação 1: Porém, se em um operação algébrica, determinada entrada estiver um precisão dupla, enquanto as demais em precisão simples, o numpy "promove" a operação para a maior precisão.
Observação 2: Os resultados apresentados foram obtidos ao executar os scripts no Python 3.8.3, numpy 1.18.5, sobre o Windows 10 em uma arquitetura x64.
Basicamente, o algoritmo desenvolvido:
Armazena \(A\) e \(b\) com precisão simples ou dupla, conforme o usuário define a variável precisão, para float32 ou float64, respectivamente. Para o exercício em questão, utilizaremos precisão simples conforme solicitado.
Concatena \(A\) e \(b\) na matriz \(G\) e realiza o processo de Eliminação Gaussiana, com pivoteamento parcial (em que apenas a coluna atual é avaliada na busca do maior pivô. Esse feature foi implementado de forma opcional, e será desabilitado em alguns ensaios mais a frente). \({G}\) é retornada então no formato triangular superior (com relação a parte quadrada referente a \({A}\))
Resolve o sistema linear a partir da última linha de \({G}\), em que é possível determinar diretamente o valor de \(x_k\), \(k\) índice da última linha, e a partir deste solucionar a linha imediatamente anterior, e assim consecutivamente.
Que é implementado por:
import numpy as np
import sys
def eliminacao_gaussiana(A, b, precisao, pivoteamento_parcial):
''' Dados A e b, realiza processo de Eliminação Gaussiana,
com ou sem pivoteamento parcial
'''
G = np.concatenate((A, b), axis=1)
n = G.shape[0]
# inicializa array de multiplicadores para as operações lineares da
# Eliminação Gaussiana:
m = np.ones(n, dtype=precisao)
for k in range(n):
if(pivoteamento_parcial == "on"):
# obtem índice da linha que contém elemento da maior valor para a
# coluna em questão.
# note que consideramos só as linhas abaixo das já processadas:
indice_linha_pivo = abs(G[k:, k]).argmax() + k
if(indice_linha_pivo != k):
linha_pivo = G[indice_linha_pivo, :].copy()
G[indice_linha_pivo] = G[k]
G[k] = linha_pivo
if(G[k, k] != 0):
for i in range(k + 1, n):
m[i] = np.divide(G[i, k], G[k, k], dtype=precisao)
for j in range(k, n + 1):
for i in range(k + 1, n):
G[i, j] = np.subtract(
G[i, j], np.multiply(m[i], G[k, j], dtype=precisao),
dtype=precisao)
else:
sys.exit("Matriz A é singular")
return G
def resolva(A, b, precisao, pivoteamento_parcial="on"):
''' Dados A e b, resolve o Sistema Linear, com pivoteamento parcial
por default, ou sem caso definido na chamada
'''
G = eliminacao_gaussiana(A, b, precisao, pivoteamento_parcial)
n = G.shape[0]
x = np.zeros((n, 1), dtype=precisao)
for i in range(n-1, -1, -1):
x[i, 0] = np.divide(G[i, n], G[i, i], dtype=precisao)
for j in range(i):
G[j, n] = np.subtract(G[j, n], np.multiply(G[j, i], x[i, 0],
dtype=precisao),
dtype=precisao)
return x
def imprima(desc, array):
''' imprime o resultado, com todos os dígitos armazenados
conforme precisão utilizada
'''
print(desc, ":")
[print("{}{} = {:0.99g}".format(desc, i[0]+1, p))
for i, p in np.ndenumerate(array)]
#%% Main
if __name__ == '__main__':
#%% a)
# define precisão simples:
precisao = "float32"
A = np.array(
[[21.0, 67.0, 88.0, 73.0],
[76.0, 63.0, 7.0, 20.0],
[ 0.0, 85.0, 56.0, 54.0],
[19.3, 43.0, 30.2, 29.4]], dtype=precisao)
b = np.array([[141.0],
[109.0],
[218.0],
[ 93.7]], dtype=precisao)
x = resolva(A, b, precisao)
print("\nitem a)\n")
imprima("x", x)
Utilizando essa rotina com precisão simples, note que obtemos o seguinte resultado:
x1 = -1.00000035762786865234375
x2 = 2.000000476837158203125
x3 = -3.000000476837158203125
x4 = 4
E assim podemos notar a imprecisão com relação a resposta esperada (\(x=[-1, 2, -3, 4]\)).
Item b)
O residual \({r=b-Ax}\) é então calculado por meio de aritimética de precisão dupla, com o resultado sendo armazenado inicialmente na variável \({r_double}\), e então convertido em um vetor de precisão simples, \({r}\), conforme solicitado pela questão, pelo segmento de código abaixo:
#%% b)
# calcula o resíduo r = b - Ax, com precisão dupla:
r_double = np.subtract(b, np.matmul(A, x, dtype="float64"), \
dtype="float64")
# converte o resultado, r_double, para precisão simples, e armazena em r:
r = r_double.astype(precisao)
# imprime os desvios encontrados e a norma 1 de r:
print("\nitem b)\n")
imprima("r", r)
print("||r||_1 =", np.linalg.norm((r), ord=1))
Os valores obtidos são:
r1 = 1.752376556396484375e-05
r2 = 4.76837158203125e-07
r3 = -1.3828277587890625e-05
r4 = 7.9870233093970455229282379150390625e-07
Com norma 1 (para efeitos de comparação posterior) igual a:
||r||_1 = 3.2627584e-05
Item c)
Resolve-se então o sistema linear \({Az=r}\), obtendo-se a solução "aprimorada", \({x+z}\), conforme segmento de código a seguir:
#%% c)
# encontra z tal que Az = r, realizando o cálculo com precisão simples:
z = resolva(A, r, precisao)
# calcula a solução "aprimorada", x + z:
x_apr = np.add(x, z, dtype=precisao)
print("\nitem c)\n")
imprima("z", z)
print("\nApós aprimoramento:")
imprima("x", x_apr)
Em que \({z}\) obtido é:
z1 = 3.5762786865234375e-07
z2 = -4.76837158203125e-07
z3 = 4.76837158203125e-07
z4 = 0
E a solução aprimorada é:
x1 = -1
x2 = 2
x3 = -3
x4 = 4
Assim, com apenas um "aprimoramento", chega-se na solução exata por meio do procedimento descrito até então.
Para que seja possível explorar o item d da questão, ou seja, obter-se mais passos de refinamento, tentamos então reduzir a precisão do processo de Eliminação Gaussiana, "desligando" o passo de pivoteamento parcial, e assim mostrando que dessa forma realmente pioramos a precisão da resposta. Para isso, executa-se o seguinte segmento de código:
#%% c) alternativa
x = resolva(A, b, precisao, pivoteamento_parcial="off")
# realiza a operação r = b - Ax:
r = np.subtract(b, np.matmul(A, x, dtype="float64"), dtype="float64") \
.astype(precisao)
# encontra z tal que Az = r, realizando o cálculo com precisão simples:
z = resolva(A, r, precisao, pivoteamento_parcial="off")
# calcula a solução "aprimorada", x + z:
x_apr = np.add(x, z, dtype=precisao)
print("\nitem c) sem pivoteamento parcial:\n")
imprima("x", x)
imprima("r", r)
print("||r||_1 =", np.linalg.norm((r), ord=1))
imprima("z", z)
print("\nApós aprimoramento:")
imprima("x", x_apr)
Como resultado "pré aprimoramento", obtemos:
x1 = -0.99985682964324951171875
x2 = 2.0006465911865234375
x3 = -2.9977271556854248046875
x4 = 3.996625423431396484375
r1 = 5.60283660888671875e-06
r2 = -3.45706939697265625e-05
r3 = -1.239776611328125e-05
r4 = 6.04099886913900263607501983642578125e-06
Podemos notar que a norma 1 do resíduo é quase o dobro da solução que utiliza pivotamento parcial:
||r||_1 = 5.8612295e-05
E após resolvido \({Az=r}\), a solução "aprimorada" \({x+z}\) é:
x1 = -0.999999940395355224609375
x2 = 2.0000002384185791015625
x3 = -2.99999904632568359375
x4 = 3.9999988079071044921875
Ou seja, ainda há oportunidade de melhoria.