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

标签： pythonarraysalgorithmnumpyindex-error

Answer 1

这是因为维基百科文章中的定义方式a, b, c和d定义。如果你仔细看，你会发现 is 的第一个元素，afora2的循环也是newDto 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