Блог

Я работаю над реализацией алгоритма устранения Гаусса в python с использованием numpy . Работая над этим, я заметил странное поведение. Вот что я сделал до сих пор:

 def gauss_elimination(a, b):
    n = len(b)
    print(a)
    for k in range(0, n-1):
        for i in range(k 1, n):
            lam = a[i, k] / a[k, k]
            a[i, k:n] = a[i, k:n] - (lam * a[k, k:n])
            b[i] = b[i] - lam * b[k]
            print(a)
    return b
 

С помощью этого кода, учитывая следующие массивы:

 a = np.array([[4, -2, 1], [-2, 4, -2], [1, -2, 4]])
b = np.array([11, -16, 17])
 

В результате будет:

 array([ 11, -10,  10])
 

Вот как алгоритм изменил массив a:

 array([[ 4, -2,  1],
       [ 0,  3, -1],
       [ 0,  0,  2]])
 

Что неверно. По какой-то причине значение во второй строке и третьем столбце равно -1, когда должно быть -1,5. Я вставил некоторую печать, чтобы увидеть, что на самом деле происходило, и по какой-то причине numpy усекает результат. Это измененный код:

 def gauss_elimination(a, b):
    n = len(b)
    print(a)
    for k in range(0, n-1):
        for i in range(k 1, n):
            lam = a[i, k] / a[k, k]
            print(a[i, k:n])
            print(lam * a[k, k:n])
            print(a[i, k:n] - lam * a[k, k:n])
            a[i, k:n] = a[i, k:n] - (lam * a[k, k:n])
            b[i] = b[i] - lam * b[k]
            print(a)
    return b
 

И учитывая те же массивы, которые были определены некоторое время назад, результаты будут:

 [[ 4 -2  1]
 [-2  4 -2]
 [ 1 -2  4]]
[ 0.   3.  -1.5] # This shows that the value is being calculated as presumed
[[ 4 -2  1]
 [ 0  3 -1] # But when the object a is updated, the value is -1 and not -1.5
 [ 1 -2  4]]
[ 0.   -1.5   3.75]
[[ 4 -2  1]
 [ 0  3 -1]
 [ 0 -1  3]]
[0.                 2.6666666666666665]
[[ 4 -2  1]
 [ 0  3 -1]
 [ 0  0  2]]
 

Я немного сбит с толку. Возможно, я мог ошибиться, но печать показывает, что все рассчитано так, как должно быть. Какие-нибудь советы?

Проблема

Ваш массив dtype -это "int32" :

 >>> import numpy as np
>>> x = np.array([1, 2, 3])
>>> x[1] / x[2]
0.6666666666666666
>>> x[0] = x[1] / x[2]
>>> x
array([0, 2, 3])
>>> x.dtype
>>> dtype('int32')
 

«Решение»

Вы можете использовать dtype="float64" создание экземпляра at для устранения этой проблемы:

 >>> x = np.array([1, 2, 3], dtype="float64")
>>> x[0] = x[1] / x[2]
>>> x
array([0.66666667, 2.        , 3.        ])
 

Конечно, вы также можете сделать это без необходимости явно указывать dtype , создав экземпляр своего массива с плавающими символами в первую очередь:

 >>> x = np.array([1., 2., 3.])
>>> x.dtype
dtype('float64')
 

Однако всякий раз, когда вы имеете дело с численными методами программирования, вы должны знать об ограничениях арифметики с плавающей запятой. Взгляните на примечание в нижней части этого ответа для получения более подробной информации.

Объяснение

Массивы NumPy являются однородными массивами, что означает, что каждый элемент относится к одному и тому же типу. dtype Или «тип данных» определяет, сколько памяти требуется каждому элементу массива. Это, в свою очередь, определяет объем памяти всего массива в памяти.

Поскольку массивы NumPy однородны, при создании массива он хранится в непрерывном блоке памяти, что делает доступ к элементам очень обширным и легким, поскольку все элементы занимают одинаковое количество байтов. Это одна из причин, по которой NumPy так быстр.

Просто к вашему сведению, это очень упрощенное объяснение того, что здесь происходит. По сути, массивы NumPy-это массивы одного типа ( "int8" , или "int32" , или "float64" и т. Д.).

Такая операция, как деление на целые числа, всегда приводит к плавающему значению, и, поскольку NumPy не хочет предполагать, что вы хотели изменить dtype весь массив (что потребовало бы создания всего нового массива в памяти), он спокойно поддерживает dtype существующий массив.

Примечание

Всякий раз, когда вы имеете дело с численными вычислениями, вам необходимо знать об ограничениях арифметики с плавающей запятой. Арифметика с плавающей запятой подвержена неточностям из-за того, как эти числа хранятся в памяти.

В случае решения систем уравнений с использованием линейной алгебры, даже после решения для определенного вектора x с Ax = b использованием прямого алгоритма исключения Гаусса, основанного на арифметике с плавающей запятой, выражение Ax не обязательно будет точно равно b . Я рекомендую читать дальше: https://en.wikipedia.org/wiki/Floating-point_arithmetic

Кстати, вы наткнулись на дисциплину, известную как численный анализ! Более конкретно, численная линейная алгебра. Это глубокая тема. Исключение Гаусса-прекрасный алгоритм, но я также надеюсь, что этот ответ побудит будущих читателей изучить более продвинутые алгоритмы решения систем уравнений.