如何解决使用 Matlab 的反斜杠运算符来反转稀疏矩阵会导致某些条目向下舍入为零
我正在使用 Matlab 的 backslash 运算符来求解写为两个矩阵 M1
和 M2
的方程组。这两个矩阵是方阵和 tridiagonal,因此我将它们定义为稀疏矩阵。例如,每个条目的尺寸为 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,N);
M_out = M1\M2;
例如,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 中是否有更好的方法来解决稀疏矩阵求逆的问题?还是我以错误的方式使用了反斜杠运算符?
解决方法
如果矩阵的大小不会增长太多,我建议进行完整的符号计算:
N = 5;
syms a
M1 = diag(-a*ones(N-1,1),-1) + diag((1 + 2*a)*ones(N,0) + diag(-a*ones(N-1,+1);
M2 = diag(+a*ones(N-1,-1) + diag((1 - 2*a)*ones(N,0) + diag(+a*ones(N-1,+1);
M_out = M1\M2;
M_num_1e10 = subs(M_out,a,1e10);
M_num_1e16 = subs(M_out,1e16);
vpa(M_num_1e10)
vpa(M_num_1e16)
在这种情况下,您将需要 Symbolic Math Toolbox。如果您没有它,我认为您应该考虑迁移到 Python 并使用 SymPy。
编辑:
考虑到您定义问题的方式,您的计算需要更高的精度。双精度是不够的。例如,在双精度中 (1e16+1) 必须四舍五入为 (1e16),也就是说 (1e16+1)-(1e16) 等于零。所以你的问题始于矩阵的主对角线。 MATLAB 仅通过其符号工具箱提供扩展精度。
如果您想坚持使用双精度,您可以依靠所谓的双双算术自行扩展双精度。我说你必须自己做,因为我认为MATLAB没有开源的double-double库。
版权声明:本文内容由互联网用户自发贡献,该文观点与技术仅代表作者本人。本站仅提供信息存储空间服务,不拥有所有权,不承担相关法律责任。如发现本站有涉嫌侵权/违法违规的内容, 请发送邮件至 dio@foxmail.com 举报,一经查实,本站将立刻删除。