Я хочу использовать метод scipy.linalg.solve_banded для решения СЛАУ, однако не могу понять в каком виде представлена ленточная матрица для этого метода.
В моей программе ленточная матрица записывается следующим образом:
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()
И в итоге получается к примеру такая матрица
[[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]]
Метод не умеет работать с такой матрицей. Я нашел эту ссылку, привел матрицу к данному виду с помощью следующего кода: (Я создаю новую матрицу B, которая является той же матрицей A, только приведенная к такому же виду, что и по ссылке)
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
Тогда матрица А принимает более компактный вид:
[[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. ]]
Но и с этой матрицей метод не хочет вычислять СЛАУ. К какому виду мне нужно привести матрицу А, чтобы заработал метод scipy.linalg.solve_banded? Как его правильно использовать?
Если это необходимо, то вот моя реализация метода прогонки:
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)