使用本征求解线性方程组

Solving Systems of Linear Equations using Eigen

本文关键字：线性方程组      更新时间：2023-10-16

我目前正在C++进行流体模拟，算法的一部分是求解稀疏线性方程组。为此，人们建议使用Eigen库。我决定使用我编写的这个简短程序来测试它：

#include <Eigen/SparseCholesky>
#include <vector>
#include <iostream>
int main() {
std::vector<Eigen::Triplet<double>> triplets;
triplets.push_back(Eigen::Triplet<double>(0, 0, 1));
triplets.push_back(Eigen::Triplet<double>(0, 1, -2));
triplets.push_back(Eigen::Triplet<double>(1, 0, 3));
triplets.push_back(Eigen::Triplet<double>(1, 1, -2));
Eigen::SparseMatrix<double> A(2, 2);
A.setFromTriplets(triplets.begin(), triplets.end());
Eigen::VectorXd b(2);
b[0] = -2;
b[1] = 2;
Eigen::SimplicialCholesky<Eigen::SparseMatrix<double>> chol(A);
Eigen::VectorXd x = chol.solve(b);
std::cout << x[0] << ' ' << x[1] << std::endl;
system("pause");
}

它给出了这两个等式：

x - 2y = -2

3x - 2y = 2

正确的解决方案是：

x = 2

y = 2

但问题是，当程序运行时，它会输出： 0.181818 -0.727273

这是完全错误的！我已经调试了几个小时，但这是一个非常短的程序，我完全按照 Eigen 网站上的教程进行操作。有谁知道导致此问题的原因？

附言我知道我使用的类是用于稀疏矩阵的，但这些类与普通矩阵类之间的唯一区别是元素的存储方式。

SimplicialCholesky是对称正定(SPD(矩阵，你组装的矩阵甚至不是对称的。默认情况下，它只读取下三角形部分中的条目，而忽略其他条目，因此它解决了：

x + 3y = -2
3x -2y = 2

如您所注意到的，对于非对称平方问题，您需要在迭代求解器世界中使用基于 LU 或BICGSTAB的直接求解器。这在文档中进行了总结。

您应该使用能够处理非对称稀疏矩阵的求解器。另一种可能的方法是寻求不是原始系统 [A]x=b 的解决方案，而是 [A]T*[A]x=[A]T*b，其中 [A]T 代表 [A] 转置。后一个系统的矩阵是对称的和正定的(只要[A]是非奇异的(。唯一的缺点是，如果原始[A]在这个意义上不是"好"的，[A]T[A]可能相当不适应。只是旨在解决此类问题的软件示例： http://members.ozemail.com.au/~comecau/CMA_LS_Sparse.htm

相关文章：