使用 Matlab 的反斜杠运算符反转稀疏矩阵会导致某些条目向下舍入为零

Question

我正在使用Matlab的反斜杠运算符来求解写为两个矩阵M1和M2的方程组。这两个矩阵是正方形和三对角的，所以我把它们定义为稀疏矩阵。例如，每个条目的维度为5x5，它们的定义如下，每个条目中的值取决于某个常数a:

N = 5;
a =  1e10;
M1 = spdiags([-a*ones(N,1)...        % Sub diagonal
(1 + 2*a)*ones(N,1)... % Main Diagonal
-a*ones(N,1)],...      % Super diagonal
-1:1,N,N);
M2 = spdiags([+a*ones(N,1)...
(1 - 2*a)*ones(N,1)...
+a*ones(N,1)],...
-1:1,N,N);
M_out = M1M2;

例如，M1的完整形式如下：

>> full(M1)
ans =
1.0e+10 *
2.0000   -1.0000         0         0         0
-1.0000    2.0000   -1.0000         0         0
0   -1.0000    2.0000   -1.0000         0
0         0   -1.0000    2.0000   -1.0000
0         0         0   -1.0000    2.0000

现在，如果我检查结果M_out中非零条目的数量，那么我可以看到它们都是非零的，这很好：

>> nnz(M_out)
ans =
25

问题是，对于常数a的较大值，我也需要这样做。然而，例如，如果改为a=1e16，则M_out的非对角线条目会自动设置为零，可能是因为它们变得太小：

>> nnz(M_out)
ans =
5

在Matlab中有更好的方法来解决稀疏矩阵的反转问题吗？或者我使用反斜杠运算符的方式不对？

Answer 1

如果矩阵的大小没有增长太多，我建议进行完整的符号计算：

N = 5;
syms a
M1 = diag(-a*ones(N-1,1),-1) + diag((1 + 2*a)*ones(N,1),0) + diag(-a*ones(N-1,1),+1);
M2 = diag(+a*ones(N-1,1),-1) + diag((1 - 2*a)*ones(N,1),0) + diag(+a*ones(N-1,1),+1);
M_out = M1M2;
M_num_1e10 = subs(M_out,a,1e10);
M_num_1e16 = subs(M_out,a,1e16);
vpa(M_num_1e10)
vpa(M_num_1e16)

在这种情况下，您将需要"符号数学工具箱"。如果您没有，我认为您应该考虑迁移到Python并使用SymPy。

编辑：

考虑到你定义问题的方式，你的计算需要更高的精度。两倍的精度是不够的。例如，在双精度中(1e16+1(必须四舍五入到(1e16(，换句话说，(1e16+1(-(1e116(等于零。所以你的问题从矩阵的主对角线开始。MATLAB仅通过其符号工具箱提供扩展的精度。

如果你想坚持双精度，你可以依靠所谓的双精度算术来扩展双精度。我说你必须自己做，因为我不认为有一个开放源码的双双库MATLAB。