Eigen求解简单稀疏微分方程组结果错误，原因何在？

我正在使用Eigen库求解稀疏微分方程组：u''(x)=1，边界条件为u(0)=u(L)=0。以下是我的实现代码：

#include <Eigen/Sparse>
#include <vector>
#include <iostream>
#include <Eigen/IterativeLinearSolvers>
// u''(x)= 1
// Boundary conditions: u(0) = u(L) = 0. 
typedef Eigen::SparseMatrix<double> SpMat; // declares a column-major 
                                           // sparse matrix type of double
typedef Eigen::Triplet<double> T;

int main(int argc, char** argv)
{ 
  int n = 3;  // number of points
  double L=1.;
  double h=L/(n-1); // deltax=h
 
  // Assembly:
  std::vector<T> coefficients; // list of non-zeros coefficients
  Eigen::VectorXd b(n);  // the right hand side-vector resulting 
                         // from the constraints

  b.setZero();
  for(int i=0; i<n; i++)
  {
      if(i==0 || i==n-1){
        b(i)=0.;
        coefficients.push_back(T(i,i,1));
      }
      else{
      coefficients.push_back(T(i,i-1,1.));
      coefficients.push_back(T(i,i,-2.));
      coefficients.push_back(T(i,i+1,1.));
      b(i) = -h * h;}
  }
 
  SpMat A(n,n);
  A.setFromTriplets(coefficients.begin(), coefficients.end());
 
  // Solving:
  Eigen::SimplicialCholesky<SpMat> chol(A); // performs a Cholesky factorization of A
  Eigen::VectorXd x = chol.solve(b);        // use the factorization to solve 

  std::cout << "A: \n" << Eigen::MatrixXd(A) << std::endl;
  std::cout << "b:\n" << b << std::endl; 
  std::cout << "x:\n" << x << std::endl; 
  return 0;
}

运行上述代码后得到如下输出：

A:  
 1  0  0  
 1 -2  1  
 0  0  1  
b:
    0
-0.25
    0
x:
-0.0833333
 0.0833333
         0

但预期结果应为：

x:
0.0
.125
0.0

即使增大n的取值，该错误仍持续存在，请问问题出在哪里？

问题根源 #

你的有限差分格式符号错误，同时矩阵类型与求解器不匹配：

1. 方程推导符号错误：
   二阶导数的中心差分近似公式为：
   $$u''(x_i) \approx \frac{u_{i-1} - 2u_i + u_{i+1}}{h^2}$$
   代入原方程$u''(x)=1$，整理后应为：
   $$u_{i-1} - 2u_i + u_{i+1} = h^2$$
   但你给b(i)赋值为-h*h，相当于求解$u_{i-1} - 2u_i + u_{i+1} = -h^2$，和原方程完全相反，直接导致结果符号错误。
2. 求解器与矩阵不兼容：
   你使用的SimplicialCholesky只适用于正定矩阵，但当前非边界行是[1,-2,1]的矩阵是负定矩阵，这会引发求解过程中的数值问题，甚至可能导致求解失败。

修正方案 #

最小修正：修正b的符号 #

仅修改右侧向量的符号即可得到正确结果：

else{
    coefficients.push_back(T(i,i-1,1.));
    coefficients.push_back(T(i,i,-2.));
    coefficients.push_back(T(i,i+1,1.));
    b(i) = h * h; // 将-h*h改为h*h
}

推荐修正：调整矩阵符号适配Cholesky分解 #

为了让矩阵变为正定矩阵，更适合Cholesky求解器，可同时调整矩阵A的符号，最终方程与原方程等价：

else{
    coefficients.push_back(T(i,i-1,-1.));
    coefficients.push_back(T(i,i,2.));
    coefficients.push_back(T(i,i+1,-1.));
    b(i) = h * h;
}

验证结果 #

修正后，当n=3时输出的x会与预期完全一致：

x:
0
0.125
0

增大n的取值后，结果也会贴合解析解$u(x)=\frac{x(L-x)}{2}$的数值近似。

内容的提问来源于stack exchange，提问作者Eder Lima de Albuquerque