Решение разреженного Ax=b в scipy

Мне нужно решить Ax=b где A является матрицей, которая представляет метод конечных разностей для PDE. Типичный размер A для двумерной задачи около (256^2) х (256^2), и она состоит из нескольких диагоналей. Следующий пример кода, как я строю A:

N = Nx*Ny  # Nx is no. of cols (size in x-direction), Ny is no. rows (size in y-direction)

# finite difference in x-direction
up1 = (0.5)*c
up1[Nx-1::Nx] = 0
down1 = (-0.5)*c
down1[::Nx] = 0
matX = diags([down1[1:], up1[:-1]], [-1,1], format='csc')

# finite difference in y-direction
up1 = (0.5)*c
down1 = (-0.5)*c
matY = diags([down1[Nx:], up1[:N-Nx]], [-Nx,Nx], format='csc')

Добавление matX а также matY вместе получается четыре диагонали. Выше для дискретизации второго порядка. Для дискретизации четвертого порядка у меня есть восемь диагоналей. Если у меня есть вторая производная, то главная диагональ также отлична от нуля.

Я использую следующий код для фактического решения:

# Initialize A_fixed, B_fixed

if const is True: # the potential term V(x) is time-independent
    A = A_fixed + sp.sparse.diags(V_func(x))
    B = B_fixed + sp.sparse.diags(V_func(x))
    A_factored = sp.sparse.linalg.factorized(A)

for idx, t in enumerate(t_steps[1:],1):
    # Solve Ax=b=Bu

    if const in False: #
        A = A_fixed + sp.sparse.diags(V_func(x,t))
        B = B_fixed + sp.sparse.diags(V_func(x,t))

    psi_n = B.dot(psi_old)

    if const is True:
        psi_new = A_factored(psi_n)
    else:
        psi_new = sp.sparse.linalg.spsolve(A,psi_n,use_umfpack=False)

    psi_old = psi_new.copy()

У меня есть пара вопросов:

  1. Какой лучший способ решить Ax=b в scipy? Я использую spsolve в sp.sparse.linalg библиотека, которая использует LU-декомпозицию. Я пытался использовать стандарт sp.linalg.solve для плотной матрицы, но это значительно медленнее. Я также попытался использовать все другие итерационные решатели в sp.sparse.linalg библиотека, но они также медленнее (для 1000x1000 все они занимают пару секунд, по сравнению с менее чем полсекунды для spsolve, и мой A скорее всего будет намного больше). Есть ли альтернативные способы сделать вычисления?

Изменить: проблема, которую я пытаюсь решить, на самом деле зависит от времени уравнения Шредингера. Если потенциальный член не зависит от времени, то, как и предполагалось, я могу предварительно разложить матрицу A первый, чтобы ускорить код, но это не работает, если потенциал изменяется во времени, так как мне нужно изменить диагональный член обеих матриц A а также B на каждом временном шаге. Для этой конкретной проблемы, есть ли способы ускорить код, используя метод, аналогичный предварительной факторизации или другие способы?

  1. Я установил umfpack, Я пробовал настройку use_umfpack Правда и Ложь, чтобы проверить это, но на удивление use_umfpack=True занимает почти вдвое больше времени, чем use_umfpack=False, Я думал, что наличие этого пакета должно ускорить. Есть идеи, в чем дело? (PS: я использую Anaconda Spyder для запуска кода, если это имеет значение)

У меня есть использование cProfile профилировать мои коды, и почти все время тратится на эту линию с spsolve, Поэтому я думаю, что оставшаяся часть кода (инициализация матрицы / задачи) в значительной степени оптимизирована, и это часть для решения матриц, которую необходимо улучшить.

Благодарю.

0 ответов

Другие вопросы по тегам