Primeira vez aqui? Seja bem vindo e cheque o FAQ!
x

Como avaliar os erros na resolução computacional de um sistema linear por intermédio da eliminação Gaussiana?

0 votos
14 visitas
perguntada Nov 7 em Ciência da Computação por bonfim_tiago (6 pontos)  

Exercício 2.10 do livro "Scientific computing: an introductory Survey (2018)", Página 101

Compartilhe

1 Resposta

0 votos
respondida Nov 7 por bonfim_tiago (6 pontos)  

A referida questão pede para que seja considerado o seguinte sistema linear de equações:

A imagem será apresentada aqui.

A seguir, é pedido para que esse sistema seja solucionado com base em uma rotina de eliminação Gaussiana. Pede-se que sejam testados vários valores para a variável epsilon, incluindo-se o que se aproxima do erro numérico do computador "unit roundoff error". É pedido, então, que, para cada valor de epsilon, sejam calculadas a condição da matriz quadrada (lado esquerdo) do sistema e o erro relativo na estimação das soluções do sistema. Ato contínuo, a questão questiona quão precisa seria a resposta gerada pelo computador para esse sistema linear.

Bom, vamos dar uma breve revisada na teoria de matrizes antes de seguir com a apresentação do código desenvolvido. A resolução de um sistema por intermédio da eliminação Gaussiana é possível, essencialmente, porque essa eliminação vai calcular a matriz inversa da matriz quadrada que se encontra a esquerda na equação apresentada acima. A seguir, é possível multiplicar ambos os lados dessa equação por essa matriz inversa. No lado esquerdo, restara a matriz identidade juntamente com o vetor coluna que contem as variáveis (soluções originais do problema), enquanto do outro lado (lado direito da equação) vai restar uma multiplicação matricial da mesma matriz inversa pela "saída" original do sistema (vetor coluna).

Por óbvio, caso não seja possível de se obter essa matriz inversa, a solução falhará. Uma matriz que não apresenta inversa é dita singular e, usualmente, afirma-se que a matriz encontra-se mal condicionada. Nessa situação, não é possível empregar exatamente essa metodologia na resolução do problema. O "condition number" de que o exercício fala refere-se precisamente a esse condicionamento da matriz quadrada do sistema que deve ser invertida. Caso esse número seja muito grande, a matriz está rumando no caminho de se tornar singular e devemos observá-la com cuidado especial, tendo em vista tratar-se de uma situação potencialmente delicada e que pode implicar em erros na solução do problema.

Depois, observa-se que a Seção 2.3.4 do livro apresenta um limite superior para o erro relativo que se pode obter nos parâmetros solução do problema em função do número de condição dessa matriz, bem como no parâmetro epsilon do computador do usuário. Esse parâmetro epsilon se refere-se a precisão numérica das variáveis com as quais é possível se realizar contas no computador do usuário. Para um sistema em ponto flutuante que opere em 64 bits, essa variável encontra-se estimada em "\(" 2^{-52} \approx 2.22 \times 10^{-16} "\)" .

O limite superior de erro apresentado seria dado por:

"\(" \frac{||\hat{x} - x||}{||x||} \lessapprox cond(A) \epsilon_{mach} "\)"

De forma que a interpretação desse resultado, já antecipando alguns resultados, é que a solução computada por ser aguardada para ter log10(cond(A)) dígitos decimais de precisão relativo a precisão da própria entrada. Assim sendo, por exemplo, um número de condição superior a 10^4 implica que nenhum dos dígitos obtidos em nossa resposta estaria correto a não ser que os dados de entrada apresentassem uma precisão superior a esses quatro dígitos.

Bom, o código desenvolvido para revolver esse problema segue apresentado abaixo:

import numpy as np


def gepp(a, b):
    # Gaussian elimination with partial pivoting.
    # % input: A is an n x n nonsingular matrix
    # %        b is an n x 1 vector
    # % output: x is the solution of Ax=b.
    # % post-condition: A and b have been modified.
    n = len(a)
    if b.size != n:
        raise ValueError("Invalid argument: incompatible sizes between A & b.", b.size, n)
    # k represents the current pivot row. Since GE traverses the matrix in the upper
    # right triangle, we also use k for indicating the k-th diagonal column index.
    for k in range(n - 1):
        # Choose largest pivot element below (and including) k
        maxindex = abs(a[k:, k]).argmax() + k
        if a[maxindex, k] == 0:
            raise ValueError("Matrix is singular.")
        # Swap rows
        if maxindex != k:
            a[[k, maxindex]] = a[[maxindex, k]]
            b[[k, maxindex]] = b[[maxindex, k]]
        for row in range(k + 1, n):
            multiplier = a[row][k] / a[k][k]
            # the only one in this column since the rest are zero
            a[row][k] = multiplier
            for col in range(k + 1, n):
                a[row][col] = a[row][col] - multiplier * a[k][col]
            # Equation solution column
            b[row] = b[row] - multiplier * b[k]
    print("Input left side of the system:", a)
    print("Input right side of the system:", b)
    x = np.zeros(n)
    k = n - 1
    x[k] = b[k] / a[k, k]
    while k >= 0:
        x[k] = (b[k] - np.dot(a[k, k + 1:], x[k + 1:])) / a[k, k]
        k = k - 1
    return x


if __name__ == "__main__":

    # Control parameter to create the matrices
    epsilon = 1.49 * 10 ** -13

    # Taking a 2 * 2 matrix (left side of the equation)
    a_matrix = np.array([[1, 1 + epsilon],
                         [1 - epsilon, 1]])

    # Taking a 2 * 1 matrix (right side of the equation)
    b_matrix = np.array([[1 + (1 + epsilon) * epsilon], [1]])

    system_solution = gepp(np.copy(a_matrix), np.copy(b_matrix))

    print(system_solution)

    # Evaluate and print Matrix Condition Number
    a_cond_number = np.linalg.cond(a_matrix, p=None)

    print("Condition of the square matrix: ", a_cond_number)

    # Evaluate the determinant of the square matrix to check if its becoming "singular"

    a_det = np.linalg.det(a_matrix)

    print("Determinant of the square matrix: ", a_det)

    # Evaluate and print coefficients relative error in the solution
    error_var_one = np.abs(100 * (1 - system_solution[0]) / 1)
    error_var_two = np.abs(100 * (epsilon - system_solution[1]) / epsilon)

    print("Relative error in output variables: [%]", error_var_one, "and", error_var_two)

Vamos rodar essa simulação para alguns valores de epsilon

epsilon = 1

Com resultado do processamento:

Input left side of the system: [[1 2]
[0 1]]
Input right side of the system: [[3]
1]
[1. 1.]
Condition of the square matrix: 5.828427124746189
Determinant of the square matrix: 1.0
Relative error in output variables: [%] 0.0 and 0.0

Nessa situação, não houve erro nenhum na estimativa da solução. Além disso, como o log10(5.8284) = 0,76, de forma que bastaria ao vetor de entrada conter um único dígito de precisão numérica para que todos os dígitos da solução fossem confiáveis.

epsilon = 1000

Com resultado do processamento:

Input left side of the system: [[-999 1]
[ 0 1001]]
Input right side of the system: [[ 1]
[1001001]]
[ 1.000001 1000.000999]
Condition of the square matrix: 1.0020020010000001
Determinant of the square matrix: 999999.9999999995
Relative error in output variables: [%] 0.0001000001000006634 and 9.99000998945121e-05

Nessa situação, houve um pequeno erro na solução do sistema, mas que se encontra dentro do limite apresentado pela equação de erro do livro.

Bom, mas se não temos erros significativos para valores "suficientemente" grandes de epsilon, vamos reduzir esses valores agora.

epsilon = 0.00001

Com resultado do processamento:

Input left side of the system: [[1.00000000e+00 1.00001000e+00]
[9.99990000e-01 9.99998973e-11]]
Input right side of the system: [[1.0000100e+00]
[8.8817842e-16]]
[1.00000112e+00 8.88179332e-06]
Condition of the square matrix: 40000084898.939415
Determinant of the square matrix: 9.99998972517348e-11
Relative error in output variables: [%] 0.00011182178591973013 and 11.182066771101763

Nessa situação, o erro na estimativa aumentou significativamente. Como pode ser observado, a condição da matriz explodiu e o seu determinante está começando a se aproximar de zero, o que tornaria a matriz singular.

Conforma orientado pelo livro, vamos pegar na próxima etapa de simulação um valor de "\("\epsilon = \sqrt{\epsilon_{mach}} \approx 1.49 \times 10^{ -13}"\)"

E os resultados da simulação seguem:

Input left side of the system: [[1. 1.]
[1. 0.]]
Input right side of the system: [[1.]
[0.]]
[nan nan]
Condition of the square matrix: inf
Determinant of the square matrix: 0.0
Relative error in output variables: [%] nan and nan

Em que vamos observar que a matriz se tornou singular, como consequência, seu número de condição agora tornou-se ilimitado e não mais é possível solucionar esse problema linear por intermédio da referida metodologia, de sorte que as estimativas para as soluções do sistema são, agora, nan = not a number.

Todos os resultados obtidos são coerentes com aqueles apresentados no referido livro.

comentou Nov 11 por Renan Abrantes (16 pontos)  
Thiago,

Excelente resolução! Você especificou muito bem a questão, passando por todos os pontos perguntados pelos autores do livro. Além disso, contextualizou o significado de eliminação gaussiana, explicando ponto a ponto o que é feito para se chegar aos resultados.
Quanto ao código escrito, achei ele bastante claro e organizado. Consegui entender de primeira o que é feito. Os resultados também saem de forma muita explícita, conforme o que se espera do exercício. Conseguir reproduzir os resultados. de maneira idêntica ao que você escreveu aqui.
Como eu não tenho absolutamente nada a adicionar à sua resolução, farei apenas um comentário lúdico: a primeira resolução pelo método de eliminação a que se tem acesso é proveniente de um livro chinês, "Os Nove Capítulos da Arte Matemática", cuja última "edição" é do século 2 D.C. Esse livro se baseia em métodos mais gerais de se resolver um problema, em contraponto à abordagem axiomática dos matemáticos gregos. Colocando em outra perspectiva, podemos afirmar que a eliminação gaussiana na verdade pode ter nascido muito antes de Gauss, originando-se de uma abordagem indutiva em vez de uma abordagem dedutiva!
...