0

Я хочу использовать метод 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)

1 ответ 1

1

Проблема в том как ты составляешь передаваемую матрицу. Она у тебя по размерностям не подходит. Ссылку не смотрел, но в документации написано как нужно составлять передаваемую матрицу. https://docs.scipy.org/doc/scipy/reference/generated/scipy.linalg.solve_banded.html

ab[u + i - j, j] == a[i,j]

Вот по этой формуле. Нужно пробежаться двумя циклами и составить эту новую матрицу с размерами (l+u+1,M).

Там еще есть звездочки. Я так понимаю это не инициализированные элементы матрицы ab.

M - количество столбцов(для индекса j).

u (upper) - параметр количество не нулевых диагоналей НАД главной диагональю.

l (lower) - параметр количество не нулевых диагоналей ПОД главной диагональю.

Пример там есть, но там нет этих циклов по перестройке матрицы "a" в матрицу "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]

Ваш ответ

Нажимая на кнопку «Отправить ответ», вы соглашаетесь с нашими пользовательским соглашением, политикой конфиденциальности и политикой о куки

Всё ещё ищете ответ? Посмотрите другие вопросы с метками или задайте свой вопрос.