简单稀疏方程组的特征解是错误的。为什么?

问题描述 投票:0回答: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 的大小,错误仍然存在。

c++ sparse-matrix eigen finite-difference
© www.soinside.com 2019 - 2024. All rights reserved.