Решение системы линейных уравнений с использованием spsolve - PullRequest
Новая
       30

Решение системы линейных уравнений с использованием spsolve

0 голосов
/ 16 декабря 2018

У меня очень большая система уравнений, согласно которой диагонали матрицы +1, +2, -1, -2 заполнены, а все остальные элементы равны нулю (кроме трех элементов в верхнем правом и нижнем левом углу).угол матрицы).Я попытался с помощью a-решить, но это не дает правильный результат.Могу ли я получить помощь относительно того, как правильно использовать это, или если есть более эффективный алгоритм для этого типа разреженной матрицы? 

from scipy.sparse import csc_matrix
from scipy.sparse.linalg import spsolve

def single_solition(x, alpha, shift, dx):
    return 12*((alpha*(1/np.cosh(alpha*(x-shift))))**2)
    
dx = 0.1
extent = 20
x = np.arange(-1*extent, extent, dx, dtype = float)
z = single_solition(x, 1, 0, dx)

def construct_Jacobian(dim, dx,z):
    J = np.zeros((dim, dim))
    
    minus_one = (z/(2*dx)) - (1/(dx**3))

    np.fill_diagonal(J[1:], minus_one)
    plus_one = -1*(z/(2*dx)) + (1/(dx**3))

    np.fill_diagonal(J[:,1:], plus_one[1:])
    
    minus_two = 0.5*(dx**3)
    np.fill_diagonal(J[2:], minus_two)
    np.fill_diagonal(J[:,2:], -minus_two)
    
    J[0, -2], J[1, -1] = minus_two, minus_two
    J[-1, 1], J[-2, 0] = -minus_two, -minus_two
    
    J[0, -1] = minus_one[-1]
    J[-1, 0] = plus_one[0]

    return J

J = construct_Jacobian(len(z), dx, z)
A = csc_matrix(J, dtype=float)
B = csc_matrix(np.transpose([-1*z]), dtype=float)
x = spsolve(A, B)
print(A.todense()*np.transpose(np.matrix(x)) + z)
...