I can't figure out why this code is not working:
#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;
}
When I run it I get the following output:
A:
1 0 0
1 -2 1
0 0 1
b:
0
-0.25
0
x:
-0.0833333
0.0833333
0
What is wrong? The result should be
x:
0.0
.125
0.0
If I increase the size of n, the error persists.