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, 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

Answer 2

向量 a 和 c 的长度应与 b 和 d 的长度相同，因此只需 prepend/append 分别将它们置零即可。此外，范围应该是 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)

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

what is the solution to a tridiagonal matrix in python?

python

arrays

algorithm

numpy

index-error