稀疏对角矩阵求解器

Question

我希望在MatLab中解决一个线性系统（对应于一个写成有限差分格式的PDE系统的两个方程）。该系统矩阵的操作（对应于PDE系统的一种扩散项）可表示为符号形式（其中u是一个未知场，n是时间步长，j是网格点）： 完整表达式如下： 上述矩阵应被视为A，其中A * U^n+1 = B 是该系统。 U交替包含'u'和'v'（PDE系统的第二个未知场）：U=[u_1,v_1,u_2,v_2,...,u_J,v_J]。 到目前为止，我一直在以以下昂贵的方式使用spdiags和diag填充该矩阵：

    E=zeros(2*J,1);

    E(1:2:2*J) = 1;
    E(2:2:2*J) = 0;

    Dvec=zeros(2*J,1);

        for i=3:2:2*J-3
                 Dvec(i)=D_11((i+1)/2);    
        end

        for i=4:2:2*J-2
                 Dvec(i)=D_21(i/2);
        end

    A = diag(Dvec)*spdiags([-E,-E,2*E,2*E,-E,-E],[-3,-2,-1,0,1,2],2*J,2*J)/(dx^2);`

而对于解决方案：

[L,U]=lu(A);
 y = L\B; 
 U(:) =U\y; 

其中B是右手边的向量。

显然，这种方法非常耗费资源，因为需要建立一个JxJ的矩阵，进行JxJ矩阵乘法等操作。

那么我的问题来了：是否有一种方法可以在不传递矩阵的情况下解决该系统，例如通过传递向量Dvec或直接传递D_11和D_22？这将节省大量内存和处理时间！

- usumdelphini

Answer 1

Matlab不会将稀疏矩阵存储为JxJ数组，而是作为大小为O(J)的列表。详见http://au.mathworks.com/help/matlab/math/constructing-sparse-matrices.html。由于您正在使用spdiags函数构造A，Matlab应该已经识别A为稀疏矩阵，如果在控制台视图中显示A，您应该确实看到这样的列表。

对于像您的三对角矩阵这样的矩阵，L和U矩阵应该已经是稀疏的。

因此，您只需要确保\操作符根据http://au.mathworks.com/help/matlab/ref/mldivide.html中的规则使用相应的稀疏算法。不清楚向量B是否已被认为是稀疏的，但您可以将其重新转换为对角矩阵，这应该肯定被认为是稀疏的。

- Dave Kielpinski