(IndexError) n x n 矩阵的反向替换

Question

我尝试使用 Python 实现反向替换，从 Lx = b 解决 x 但出现错误。

而且我不知道如何进行迭代以了解问题所在。

我是不是误解了一个概念？

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]