(IndexError) Обратная замена матрицы nxn

Question

Я пытаюсь реализовать Back-Substitution с помощью Python, решая x из Lx = b , но я получил ошибку.

И я не знаю, как ямогу делать итерации, чтобы понять, что не так.

Не понимаю ли я концепцию?

import numpy as np
def backSub(L: np.ndarray, b: np.ndarray)-> np.ndarray:
    length = len(L)
    x = np.zeros(length)
    for k in range(length - 1, -1, -1):
        x[k] = (b[k] - np.dot(L[k, k + 1:], x[k + 1:])) / L[k, k]

    return x


M = np.array([[2, 0,0],
             [1, 3, 0],
             [2, 3, 4]])

s = np.array([[2, 2, 0]])

print(backSub(M, s))

когда я запускаю свою Программу, я получаю следующее:

Вывод

x[k] = (b[k] - np.dot(L[k, k + 1:], x[k + 1:])) / L[k, k]
IndexError: index 2 is out of bounds for axis 0 with size 1

ожидаемый результат

Просто вектор x от ** L x = b **

Answer 1

Вы можете сделать исключение и обратную замену в той же функции.

import numpy as np

def gauss_elimin(a,b):
    """
    solve [a]{x} = {b} for {x}
    """
    a = np.copy(a) # copy a and b if you want them for later purpose
    b = np.copy(b) #copying will not overwrite original a and b [id(a) is different before and after copy]
    n = len(b)
    assert (np.all(np.shape(a) ==(n,n))) # check if a is a square matrix
    # Elimination Phase
    for k in range(0,n-1): # pivot row
        for i in range(k+1,n): #  rows being transformed
            if a[i,k] != 0.0: # skip if a(i,k) is already zero
                lam = a [i,k]/a[k,k]
                a[i,k+1:n] = a[i,k+1:n] - lam*a[k,k+1:n]
                b[i] = b[i] - lam*b[k]
    # Back substitution
    for k in range(n-1,-1,-1):
        b[k] = (b[k] - np.dot(a[k,k+1:n],b[k+1:n]))/a[k,k]
    return b # b now contains the solution x

Что может быть неправильным в вашей функции, так это то, что вам не хватает n

x[k] = (b[k] - np.dot(L[k, k + 1:n], x[k + 1:n])) / L[k, k]

Проверка функции

a = np.random.rand(4,4)*10
b = np.random.rand(4)*10

x1 = gauss_elimin(a, b)

from scipy.linalg import solve

x2 = solve(a, b)

print(np.allclose(x1, x2)) >> prints True

Изменение вашей функции

def back_subs(u, b):
    u = u.copy()
    b = b.copy()
    n = b.size

    x = np.zeros(n)

    for k in range(n-1,-1,-1):
        x[k] = (b[k] - np.dot(u[k,k+1:n],x[k+1:n]))/u[k,k]

    return x

Проверка вашей функции после модификации

a = np.random.rand(4, 4)
b = np.random.rand(4)

a_u = np.triu(a)

x1 = back_subs(a_u, b)
x2 = solve(a_u, b)

print(x1)
print(x2)

[ 0.55056814 -1.32860072  0.24254532  1.49865603]
[ 0.55056814 -1.32860072  0.24254532  1.49865603]