python中三对角矩阵的解决方案是什么？

Question

我现在正在处理一个关于三对角矩阵的问题，我使用了 wiki https://en.wikipedia.org/wiki/Tridiagonal_matrix_algorithm 中的三对角矩阵算法来实现一个解决方案，我已经尝试过，但我的解决方案并不完整。

我很困惑，我需要帮助，这也是我使用 jupyter notebook 的脚本

import numpy as np 

# The diagonals and the solution vector 

b = np.array((5, 6, 6, 6, 6, 6, 6, 6, 5), dtype=float) # Main Diagonal
a= np.array((1, 1, 1, 1, 1, 1, 1, 1), dtype = float) # Lower Diagonal
c = np.array((1, 1, 1, 1, 1, 1, 1, 1), dtype = float) # Upper Diagonal
d = np.array((3, 2, 1, 3, 1, 3, 1, 2, 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 -1):
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

Answer 1

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

import numpy as np 

# The diagonals and the solution vector 

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

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

#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):
    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)

Answer 2

这是因为a, b, c 和d 在维基百科文章中的定义方式。如果你仔细看，你会发现a 的第一个元素是a2，newD 的循环也是n。因此，为了让数组具有可理解的索引，我建议您添加一个幻像元素 a1。你得到：

import numpy as np 

# The diagonals and the solution vector 

b = np.array((5, 6, 6, 6, 6, 6, 6, 6, 5), dtype=float) # Main Diagonal
a = np.array((np.nan, 1, 1, 1, 1, 1, 1, 1, 1), dtype = float) # Lower Diagonal
                                                              # with added a1 element
c = np.array((1, 1, 1, 1, 1, 1, 1, 1), dtype = float) # Upper Diagonal
d = np.array((3, 2, 1, 3, 1, 3, 1, 2, 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