这是因为Wikipedia文章中如何定义a,b,c和d。如果仔细看，您会发现a的第一个元素是a2，并且newD的循环也转到n。因此，为使数组具有可理解的索引，我建议您添加一个幻像元素a1。然后您得到：

import numpy as np 

# The diagonals and the solution vector 

b = np.array((5,6,5),dtype=float) # Main Diagonal
a = np.array((np.nan,1,1),dtype = float) # Lower Diagonal
                                                              # with added a1 element
c = np.array((1,dtype = float) # Upper Diagonal
d = np.array((3,2,3,3),dtype = float) # Solution Vector

#number of rows 
n = d.size

newC = np.zeros(n,dtype= float)

newD = np.zeros(n,dtype = float)

x = np.zeros(n,dtype = float)

# the first coefficents 
newC[0] = c[0] / b[0]
newD[0] = d[0] / b[0]

for i in range(1,n - 1):
    newC[i] = c[i] / (b[i] - a[i] * newC[i - 1])

for i in range(1,n):  # Add the last iteration `n`
    newD[i] = (d[i] - a[i] * newD[i - 1]) / (b[i] - a[i] * newC[i - 1])

x[n - 1] = newD[n - 1]

for i in reversed(range(0,n - 1)):
    x[i] = newD[i] - newC[i] * x[i + 1]


x

向量a和c的长度应与b和d的长度相同，因此只需在它们的前面分别添加零即可。此外，范围应为range(1,n)，否则，您的最后一个解决方案元素应为0。您可以看到此修改后的代码，以及与已知算法的比较，表明它得到了相同的答案。

import numpy as np 

# The diagonals and the solution vector 

b = np.array((5,dtype=float) # Main Diagonal
a= np.array((0,dtype = float) # Lower Diagonal
c = np.array((1,0),dtype = float) # Solution Vector

print(b.size)
print(a.size)

#number of rows 
n = d.size

newC = np.zeros(n,n):
    newC[i] = c[i] / (b[i] - a[i] * newC[i - 1])

for i in range(1,n):
    newD[i] = (d[i] - a[i] * newD[i - 1]) / (b[i] - a[i] * newC[i - 1])

x[n - 1] = newD[n - 1]

for i in reversed(range(0,n - 1)):
    x[i] = newD[i] - newC[i] * x[i + 1]

# Test with know algorithme
mat = np.zeros((n,n))
for i in range(n):
    mat[i,i] = b[i]
    if i < n-1:
        mat[i,i+1] = a[i+1]
        mat[i+1,i] = c[i]

print(mat)

sol = np.linalg.solve(mat,d)
print(x)
print(sol)