Описание метода¶

Даны матрица $A$ размером $N\times N$, вектор $B$ высотой $N$. Требуется найти вектор $X$ такой, что $AX=B$.

Без терминологии линейной алгебры это записывается, как:

$ \begin{array}{lcr} a_{11}x_1+\ldots+a_{1N}x_N &=& b_1 \\ \ldots & & \\ a_{N1}x_1+\ldots+a_{NN}x_N &=& b_N \\ \end{array} $

Прямой ход¶

Из второго уравнения вычитаем первое, умноженное на $\frac{a_{21}}{a_{11}}$. Из третьего вычитаем первое, умноженное на $\frac{a_{31}}{a_{11}}$. И т.д. до $N$-го. Получаем систему уравнений следующего вида.

$ \begin{array}{lcr} a_{11}&x_1 +& a_{12} x_2 +\ldots+a_{1N}x_N &=& b_1 \\ 0 &x_1 +& a'_{22} x_2 +\ldots+a'_{1N}x_N &=& b'_1 \\ & \ldots & & \\ 0 &x_1 +& a'_{N2} x_2 +\ldots+a'_{NN}x_N &=& b'_N \\ \end{array} $

Затем из третьего уравнения вычитаем второе, умноженное на $\frac{a_{32}}{a_{22}}$.

И т.д., пока не выполним такое же действие над последней строкой.

В итоге мы получим систему уравнений следующего вида.

$ \begin{array}{lcr} a_{11}& x_1 +& a_{12} & x_2 + & a''_{13} & x_3 + \ldots + a_{1N} x_N &=& b_1 \\ 0 & x_1 +& a''_{22} & x_2 + & a''_{23} & x_3 + \ldots + a''_{2N} x_N &=& b''_1 \\ 0 & x_1 +& 0 & x_2 + & a''_{33} & x_3 +\ldots + a''_{3N} x_N &=& b''_3 \\ 0 & x_1 +& 0 & x_2 + & 0 & x_3 + \ldots + a''_{4N} x_N &=& b''_4 \\ & & & \ldots & & & & \\ 0 & x_1 +& & \ldots & + & 0 x_{N-1} + a''_{NN} x_N &=& b''_N \\ \end{array} $

Обратный ход¶

Очевидно, $x_N = \frac{b''_N}{a''_{NN}}$.

Далее, $x_i = \frac{b''_i - \sum_{j=i+1}^{N} a''_{ij}x_j}{a_{ii}}$

Прямой и обратный ход в вектроной форме¶

Интересно? Тогда см. ниже. Неинтересно? Тоже см. ниже.

Задание¶

Предлагается реализовать следующие недостающие функции

In [4]:
import numpy
from numpy import array
from numpy.linalg import norm
from numpy.linalg import solve as solve_out_of_the_box


def gauss(a, b):
    a = a.copy()
    b = b.copy()
    N = len(b)

    def forward():
        for i in range(N - 1):
            for j in range(i + 1, N):
                a[j] -= a[i] * (a[j][i] / a[i][i])
                b[j] -= b[i] * (a[j][i] / a[i][i])

    def backward():
        x = numpy.zeros(N, dtype=float)
        x[-1] = b[-1] / a[-1][-1]
        for i in range(N - 2, 0 + 1, -1):
            s = 0
            for j in range(i + 1, N):
                s = a[i][j] * x[j]
            x[i] = (b[i] - s) / a[i][i]
        return x

    forward()
    x = backward()
    print(a)
    print(b)
    return x


a = array([
    [1.5, 2.0, 1.5, 2.0],
    [3.0, 2.0, 4.0, 1.0],
    [1.0, 6.0, 0.0, 4],
    [2.0, 1.0, 4.0, 3]
], dtype=float)

b = array([5, 6, 7, 8], dtype=float)

print(a)
print(b)

oob_solution = solve_out_of_the_box(a, b)
solution = gauss(a, b)

print(solution)
print("Макс отклонение компоненты решения:", norm(solution - oob_solution, ord=1))

[[1.5 2.  1.5 2. ]
 [3.  2.  4.  1. ]
 [1.  6.  0.  4. ]
 [2.  1.  4.  3. ]]
[5. 6. 7. 8.]
[[ 1.5         2.          1.5         2.        ]
 [ 0.         -2.          1.         -3.        ]
 [ 0.          0.          1.33333333 -4.33333333]
 [ 0.          0.          0.          6.625     ]]
[5. 6. 7. 8.]
[0.         0.         9.1745283  1.20754717]
Макс отклонение компоненты решения: 9.853773584905657

Прямой и обратный ход в векторной форме¶

Очевидно, что на прямом ходе мы оперируем целиком уравнениями (включая свободные члены): умножаем всё уравнение на число, вычитаем одно уравнение из другого и т.д., т.е. оперируем целыми уравнениями, (строками матрицы), а не отдельными коэффициентами. При вычислениях это приводит к избыточности, но в случае, если компьютер может выполнять векторные вычисления (т.е. обрабатывать строку матрицы, как одно число), это даёт большое ускорение.

Сразу скажем, что реализация NumPy «из коробки» так не умеет. «Убедить» её это делать, пусть и не «из коробки», можно, но это не тема данного блокнота.

Итак, пусть у нас $d$ уравнений. $\mathrm{AB}$ — матрица из $d$ строк и $d+1$ столбцов, у которой правый слобец состоит из свободных членов уравнений. Запишем алгоритм в «императивной» форме.

$ab_i$ — строка, $ab_{i,j}$ — элемент матрицы.

Прямой ход¶

$ab_1 \gets \frac{ ab_1 }{ ab_{1,1}}$ — делаем 1 в левом верхнем углу

Далее для строк c индексами $j \in [2, d]$: $ab_j \gets ab_j - ab_1 * ab_{j,1}$ — вычитаем первую строку (у неё уже 1 в начале), умноженную на первый элемент $j$-й строки, делая первый элемент равным 0.

$ab_2 \gets ab_2 / ab_{2, 2}$ — делаем $1$ «почти» в левом верхнем углу.

и т.д.

В итоге прямого хода получаем $\mathrm{AB}$ с единицами на главной диагонали $A$ и нулями ниже её.

Обратный ход¶

Для $i=d-1, \ldots, 1$ повторяем следующее.

$ab_i \gets ab_i - \sum_{j=i+1}^d ab_j * ab_{i, j}$.

Результат¶

По итогам обратного хода получем $\mathrm{AB}$ такую, что $A$ — единичная матрица, а $B$ — решение системы.

In [5]:
from numpy import array, concatenate
from numpy.linalg import norm
from numpy.linalg import solve as solve_out_of_the_box
from numpy import zeros

a = array([
    [1.5, 2.0, 1.5, 2.0],
    [3.0, 2.0, 4.0, 1.0],
    [1.0, 6.0, 0.0, 4],
    [2.0, 1.0, 4.0, 3]
], dtype=float)

b = array([5, 6, 7, 8], dtype=float)


def vector_gauss(a, b):
    ab = concatenate((a, array([b]).T), axis=1)  # concatenate заодно и скопирует
    d = len(ab)  # размер по старшему измерению

    print(ab)

    # прямой
    for i in range(d):
        # ab_ii = ab[i, i]
        ab[i] /= ab[i, i]
        for j in range(i + 1, d):
            ab[j] -= ab[i] * ab[j, i]

    # обратный
    for i in range(d - 1, -1, -1):
        s = zeros(len(ab[0]), dtype=float)
        for j in range(i + 1, d):
            s += ab[j] * ab[i, j]
        ab[i] -= s

    x = ab[:, -1]  # Последний столбец
    print(ab)
    print(ab[1])
    return x


solution = vector_gauss(a, b)
oob_solution = solve_out_of_the_box(a, b)

print(solution)
print("Макс отклонение компоненты решения:", norm(solution - oob_solution, ord=1))

[[1.5 2.  1.5 2.  5. ]
 [3.  2.  4.  1.  6. ]
 [1.  6.  0.  4.  7. ]
 [2.  1.  4.  3.  8. ]]
[[ 1.          0.          0.          0.          0.8490566 ]
 [-0.          1.          0.          0.          0.05660377]
 [ 0.          0.          1.          0.          0.47169811]
 [ 0.          0.          0.          1.          1.45283019]]
[-0.          1.          0.          0.          0.05660377]
[0.8490566  0.05660377 0.47169811 1.45283019]
Макс отклонение компоненты решения: 3.4139358007223564e-15

[0.8490566 0.05660377 0.47169811 1.45283019] Макс отклонение компоненты решения: 3.635980405647388e-15

Проверим «серьёзн

Подготовимся

In [3]:
from numpy.linalg import norm, det
from numpy.linalg import solve as solve_out_of_the_box
from numpy.random import uniform

N=100
SMALL = 1e-5

def test_error(solver_function):
    # Сгенерируем хорошо обусловленную систему
    while True:
        a = uniform(0.0, 1.0, (N, N))
        b = uniform(0.0, 1.0, (N,  ))

        d = det(a)
        if abs(d) <= SMALL:  # Определитель маленький — плохо
            print(f"det: {d}")
            continue  # Попробуем ещё
        if d < 0.0:  # Отрицательный — поменяем знак
            print(f"det: {d}")
            a = -a

        try:
            oob_solution = solve_out_of_the_box(a, b)
        except Exception as e:  # Всё-таки система плохая
            print(f"exc: {e}")
            continue  # Попробуем ещё

        sb = a @ oob_solution
        if norm(sb - b, ord=1) > SMALL:
            print("Bad solution...")
            continue  # Всё же не очень система получилась =)
        
        break # Всё, считаем, что мы случайно сгенерировали хорошую систему

    tested_solution = solver_function(a, b)
    return norm(tested_solution - oob_solution, ord=1)

И собственно проверим

In [4]:
test_error(vector_gauss)

det: -2.674273033732956e+25
[[-0.15901379 -0.08900197 -0.49032723 ... -0.76554311 -0.57844704
   0.07984708]
 [-0.51041006 -0.81906461 -0.3153571  ... -0.00204936 -0.91956286
   0.9099387 ]
 [-0.13970478 -0.54434794 -0.16491843 ... -0.20210984 -0.29202303
   0.73945634]
 ...
 [-0.29665502 -0.64123577 -0.14988566 ... -0.14528975 -0.00560326
   0.53633483]
 [-0.30411977 -0.61095607 -0.12602411 ... -0.13559222 -0.5558862
   0.21794898]
 [-0.19694385 -0.53743203 -0.3498482  ... -0.86273484 -0.16035124
   0.7605858 ]]
[[  1.           0.55971227   3.08355157 ...   4.81431889   3.63771608
   -0.50213935]
 [ -0.           1.          -2.35950787 ...  -4.6031328   -1.75702275
   -1.22546679]
 [ -0.          -0.           1.         ...   2.0086858    0.72283139
   -0.1175619 ]
 ...
 [ -0.          -0.          -0.         ...   0.          39.03762941
  -11.864914  ]
 [  0.           0.           0.         ...   1.           5.00232947
   -1.66767861]
 [  0.           0.           0.         ...   0.          -7.50225532
    2.61781578]]
[-0.          1.         -2.35950787 -2.11966581 -2.21789575 -0.85477161
 -2.06325372  1.01203486 -0.58703034 -4.64152003 -1.98357312 -4.4889145
 -3.40104908 -4.62928806 -0.30236399 -2.4298808   0.05892874 -0.30463823
 -4.1971761  -4.6712647  -1.36642345 -0.79837582 -1.06805673 -0.24450619
  0.44394906 -3.80457695 -0.64288125 -2.59465664 -2.2679029  -0.34749562
 -4.7606305  -0.94735171 -2.06808221 -2.0429092   0.53013921 -2.00708798
 -4.56239339 -0.5538983  -0.79097637  0.58279885 -4.07528518 -1.21037251
 -4.12879771 -1.30396376  0.40940847 -2.49058325 -3.04998029 -5.5520489
 -1.58405574 -1.78261429 -5.20768174 -4.44199609 -3.84828975  1.06262537
 -4.54451452 -5.16061308 -2.60986029 -4.2500466  -2.5791183  -1.80623159
 -1.19457602 -2.38780967 -3.48098183 -0.55304302 -4.70371943 -4.26199238
 -3.59481553  1.10111605 -3.32446132  0.71990755 -0.63798186 -1.1307071
 -0.39108877 -4.38227618 -0.88322765  0.01623957 -1.21387297 -5.01841405
 -0.21881961 -3.35326973 -2.80435245 -1.13990976 -0.44585155 -4.70260776
 -1.81332971 -1.0993462  -3.43288159 -0.00849481 -3.84466353 -0.93958392
 -3.07925756 -2.10872723 -0.72477821  1.40924716 -2.44844715 -2.00477945
 -1.26734438 -4.39935891 -4.6031328  -1.75702275 -1.22546679]
Out[4]:
811.2469885473931