scipy.linalg.solve banded. Which matrix should I use for the run-through method?
I want to use the scipy.linalg.solve_banded method to solve SLOUGH, but I can't figure out how the ribbon matrix is represented for this method.
In my program, the tape matrix is written as follows:
A = np.zeros((n, n))
for a in range(0,n):
for b in range(a - 1,a + 2):
if (b == -1):
A[a][0] = random.random()
continue
if (b == n):
break
A[a][b] = random.random()
And the result is, for example, such a matrix
[[0.16152984 0.63902563 0. 0. 0. 0. ]
[0.628172 0.31132042 0.52016462 0. 0. 0. ]
[0. 0.54531714 0.4615466 0.93090818 0. 0. ]
[0. 0. 0.21819913 0.23413849 0.50466658 0. ]
[0. 0. 0. 0.88416538 0.74770314 0.64462594]
[0. 0. 0. 0. 0.33679632 0.58161105]]
The method does not know how to work with such a matrix. I found this link, brought the matrix to this view using the following code: (I'm creating a new matrix B, which is the same matrix A, only reduced to the same form as by reference)
B = np.random.rand(n,3)
for i in range(0,n):
for j in range(0,3):
y = i+j-1
if y >= 0 and y < n:
B[i][j] = A[i][y]
else:
B[i][j] = 0
Then the matrix A takes a more compact form:
[[0. 0.16152984 0.63902563]
[0.628172 0.31132042 0.52016462]
[0.54531714 0.4615466 0.93090818]
[0.21819913 0.23413849 0.50466658]
[0.88416538 0.74770314 0.64462594]
[0.33679632 0.58161105 0. ]]
But even with this matrix, the method does not want to calculate SLOUGH. What form do I need to bring matrix A to in order for the scipy.linalg.solve_banded method to work? How to use it correctly?
If necessary, here is my implementation of the run method:
import numpy as np
import scipy.linalg as sl
import random
#Ввод данных
n = 6
A = np.zeros((n, n))
for a in range(0,n):
for b in range(a - 1,a + 2):
if (b == -1):
A[a][0] = random.random()
continue
if (b == n):
break
A[a][b] = random.random()
f = np.random.rand(n)
x = [0] * n
#Прямой ход
m = 1;
for i in range(1,n):
m = A[i][i - 1]/A[i-1][i-1] #m = a[i]/c[i-1];
A[i][i] = A[i][i] - m*A[i-1][i] #c[i] = c[i] - m*b[i-1]
f[i] = f[i] - m*f[i-1] #f[i] = f[i] - m*f[i-1]
#Обратный ход
x[n-1] = f[n-1]/A[n-1][n-1];
for i in range(n - 2, -1, -1):
x[i]=(f[i] - A[i][i + 1]*x[i+1]) / A[i][i]
#Вывод
print(x)
1 answers
The problem is how you compose the transmitted matrix. It doesn't fit your dimensions. I didn't look at the link, but the documentation says how to make the transmitted matrix. https://docs.scipy.org/doc/scipy/reference/generated/scipy.linalg.solve_banded.html
Ab[u + i - j, j] == a[i,j]
According to this formula. You need to run through two cycles and make this new matrix with dimensions (l+u+1, M).
There are still asterisks there. I take it it's not initialized elements of the matrix ab.
M - number of columns(for index j).
U (upper) - parameter number of non-zero diagonals ABOVE the main diagonal.
L (lower) - parameter number of non-zero diagonals UNDER the main diagonal.
The example is there, but there are no such cycles for rebuilding the matrix " a "into the matrix"ab".
from scipy import linalg
import numpy as np
if __name__ == '__main__':
# для метода прогонки. 1 над диагональю. И 1 под диагональю
u = 1
l = 1
n = 5
m = 5
a = [[2, 3, 0, 0, 0],
[1, 2, 3, 0, 0],
[0, 1, 2, 3, 0],
[0, 0, 1, 2, 3],
[0, 0, 0, 2, 3]]
b = [1, 2, 3, 4, 5]
a = np.array(a)
b = np.array(b)
ab = np.zeros((u + l + 1, m))
for j in range(m):
for i in range(n):
index = u + i - j
if 0 <= index < u + l + 1:
ab[index][j] = a[i][j]
print(ab)
ans = linalg.solve_banded((l, u), ab, b)
print(ans)
Вывод
[[0. 3. 3. 3. 3.]
[2. 2. 2. 2. 3.]
[1. 1. 1. 2. 0.]]
[-13. 9. -1. -1.33333333 2.55555556]