1

我需要解决表示 PDE 有限差分法的矩阵Ax=b在哪里。A二维问题的典型大小A约为 (256^2)x(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中解决的最佳方法是什么？我在使用 LU 分解spsolve的sp.sparse.linalg库中使用 。我尝试使用sp.linalg.solve密集矩阵的标准，但速度要慢得多。我还尝试使用sp.sparse.linalg库中的所有其他迭代求解器，但它们也较慢（对于 1000x1000，它们都需要几秒钟，而对于 不到半秒spsolve，我A的可能要大得多）。有没有其他方法可以进行计算？

编辑：我试图解决的问题实际上是时间相关的薛定谔方程。如果潜在项与时间无关，那么按照建议，我可以先对矩阵进行预分解A以加快代码速度，但是如果潜在项是随时间变化的，则这不起作用，因为我需要更改两者的对角项矩阵A和B每个时间步。对于这个特定的问题，有没有办法使用类似于预分解的方法或其他方式来加速代码？

  1. 我已经安装了umfpack。我尝试设置use_umfpackTrue 和 False 来测试它，但令人惊讶的是，它use_umfpack=True比use_umfpack=False. 我认为拥有这个软件包应该可以加快速度。知道这是怎么回事吗？（PS：我正在使用 Anaconda Spyder 运行代码，如果这有什么不同的话）

我cProfile习惯于分析我的代码，几乎所有时间都花在了spsolve. 所以我认为代码的剩余部分（矩阵/问题初始化）已经优化很多，需要改进的是矩阵求解部分。

谢谢。

4

0 回答 0