Skip to content
Snippets Groups Projects
Commit dcd1e776 authored by Reichelt, Michael's avatar Reichelt, Michael
Browse files

test with hom boundary still not working

parent c0a1b745
No related branches found
No related tags found
No related merge requests found
......@@ -52,27 +52,25 @@ Eigen::VectorXd project_p0(std::vector<double> const & t_points, std::function<d
return u;
}
Eigen::MatrixXd remove_first_column_and_row(Eigen::MatrixXd const& A){
Eigen::MatrixXd M(A.rows()-1, A.cols()-1);
for(Eigen::Index i=0; i<M.rows(); ++i){
for(Eigen::Index j=0; j<M.cols(); ++j){
M(i,j) = A(i+1,j+1);
Eigen::MatrixXd set_homogenious_boundary(Eigen::MatrixXd const& A){
Eigen::MatrixXd B = A;
for(Eigen::Index i=0; i<B.rows(); ++i){
B(i,0) = 0.;
}
for(Eigen::Index i=0; i<B.cols(); ++i){
B(0,i) = 0.;
}
return M;
return B;
}
Eigen::VectorXd remove_first_row(Eigen::VectorXd const& v){
Eigen::VectorXd u(v.size()-1);
for(Eigen::Index i=0; i<u.size(); ++i){
u(i) = v(i+1);
}
return u;
double hom_scalar_prod(Eigen::VectorXd const& v, Eigen::MatrixXd const& A, Eigen::VectorXd const& u){
Eigen::MatrixXd const B = set_homogenious_boundary(A);
return v.dot(B*u);
}
TEST_CASE("Compact perturbation: matrices", "[hilbert_transformation][compact_perturbation]") {
std::vector<double> const t_points = linspace(0., T, 10);
std::vector<double> const t_points = linspace(0., T, 100);
Eigen::VectorXd const u_vec = project_p1(t_points, u);
Eigen::VectorXd const v_vec = project_p1(t_points, v);
......@@ -81,9 +79,8 @@ TEST_CASE("Compact perturbation: matrices", "[hilbert_transformation][compact_pe
SECTION("(u, Hv)"){
double const analytic_val = -0.0328351; //from mathematica
Eigen::MatrixXd const Mat = remove_first_column_and_row(hilbert_matrices.mass);
Eigen::VectorXd const tmp = Mat* remove_first_row(u_vec);
double const numeric_val = remove_first_row(v_vec).dot(tmp);
Eigen::MatrixXd const Mat = hilbert_matrices.mass;
double const numeric_val = hom_scalar_prod(v_vec, Mat, u_vec);
REQUIRE_THAT(numeric_val, Catch::Matchers::WithinAbs(analytic_val, 1e-3));
}
SECTION("(dt u, Hv)"){
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment