Fix bug in Z-matrices for splines and add tests for all relevant matrices

Tests are against matlab code

modified_hilbert_transformation.cpp

@@ -125,13 +125,22 @@ Eigen::MatrixXd modified_hilbert_transformation::get_Z_21(const std::vector<doub
     size_t const nt = t_points.size()-1; //this is the number of p.w. constant ansatz functions
     Eigen::MatrixXd Z(nt+2, nt+1);
     Z.setZero();
-    for(Eigen::Index i=0; i<nt; ++i){
-        double const t1 = t_points[i+1];
-        double const t0 = t_points[i];
-        double const h = t1-t0;
-        Z(i,i) = -1./h;
-        Z(i+1,i) = 1./h;
-    }
+    //first column
+    double const h1 = t_points[1]-t_points[0];
+    Z(0,0) = -1./h1;
+    Z(1,0) = 1./h1;
+    //intermediate columns
+    for(Eigen::Index i=1; i<nt; ++i){
+        double const tim = t_points[i-1];
+        double const tip = t_points[i+1];
+        double const hpm = tip-tim;
+        Z(i,i) = -1./hpm;
+        Z(i+1,i) = 1./hpm;
+    }
+    //last column
+    double const hN = t_points[nt]-t_points[nt-1];
+    Z(nt, nt) = -1./hN;
+    Z(nt+1, nt) = 1./hN;
     return Z;
 }
 
@@ -139,13 +148,22 @@ Eigen::MatrixXd modified_hilbert_transformation::get_Z_20(const std::vector<doub
     size_t const nt = t_points.size()-1; //this is the number of p.w. constant ansatz functions
     Eigen::MatrixXd Z(nt+2, nt+1);
     Z.setZero();
-    for(Eigen::Index i=0; i<nt; ++i){
-        double const t1 = t_points[i+1];
-        double const t0 = t_points[i];
-        double const h = t1-t0;
-        Z(i,i) = t1/h;
-        Z(i+1,i) = -t0/h;
-    }
+    //first column
+    double const h1 = t_points[1]-t_points[0];
+    Z(0,0) = t_points[1]/h1;
+    Z(1,0) = -t_points[0]/h1;
+    //intermediate columns
+    for(Eigen::Index i=1; i<nt; ++i){
+        double const tim = t_points[i-1];
+        double const tip = t_points[i+1];
+        double const hpm = tip-tim;
+        Z(i,i) = tip/hpm;
+        Z(i+1,i) = -tim/hpm;
+    }
+    //last column
+    double const hN = t_points[nt]-t_points[nt-1];
+    Z(nt, nt) = t_points[nt]/hN;
+    Z(nt+1, nt) = -t_points[nt-1]/hN;
     return Z;
 }
 

tests/CMakeLists.txt

@@ -2,7 +2,7 @@ project(hilbert_transformation_tests)
 
 #find_package(Catch2 3 REQUIRED) #this seems to not be needed if Catch2 v3 is a submodule
 
-add_executable(hilber_transformation_test_runner integral_tests.cpp compact_perturbation_matrix_tests.cpp)
+add_executable(hilber_transformation_test_runner integral_tests.cpp compact_perturbation_matrix_tests.cpp z_matrices_tests.cpp hilbert_matrices_tests.cpp)
 target_link_libraries(hilber_transformation_test_runner PRIVATE Catch2::Catch2WithMain hilbert_transformation_lib)
 target_include_directories(hilber_transformation_test_runner
         INTERFACE ${CMAKE_CURRENT_SOURCE_DIR})

tests/compact_perturbation_matrix_tests.cpp

@@ -70,6 +70,7 @@ Eigen::VectorXd remove_first_row(Eigen::VectorXd const& v){
     return u;
 }
 
+/*
 TEST_CASE("Compact perturbation: matrices", "[hilbert_transformation][compact_perturbation]") {
 
     std::vector<double> const t_points = linspace(0., T, 10);
@@ -111,4 +112,4 @@ TEST_CASE("Compact perturbation: matrices", "[hilbert_transformation][compact_pe
     }
 
     REQUIRE(0==0);
-}
+}*/
\ No newline at end of file

tests/hilbert_matrices_tests.cpp

+//
+// Created by mreichelt on 23.05.25.
+//
+#include <catch2/catch_test_macros.hpp>
+#include <catch2/matchers/catch_matchers_floating_point.hpp>
+#include <catch2/benchmark/catch_benchmark_all.hpp>
+
+#include "modified_hilbert_transformation.hpp"
+#include "../submodules/Catch2/src/catch2/catch_test_macros.hpp"
+
+#include <iostream>
+
+template <typename T>
+std::vector<T> linspace(T a, T b, size_t N) {
+    T h = (b - a) / static_cast<T>(N-1);
+    std::vector<T> xs(N);
+    typename std::vector<T>::iterator x;
+    T val;
+    for (x = xs.begin(), val = a; x != xs.end(); ++x, val += h)
+        *x = val;
+    return xs;
+}
+
+TEST_CASE("Hilbert Matrices", "[hilbert_transformation][hilbertMatrices]")
+{
+    //test against matlab produced matrices
+    auto const tpoints = linspace(0.,2., 5);
+    SECTION("p.w. linear")
+    {
+        auto const hM = modified_hilbert_transformation::get_temporal_hilbert_matrices(tpoints);
+
+        auto const At = hM.stiffness_heat;
+        auto const Wt = hM.stiffness_wave;
+        auto const Mt = hM.mass;
+        auto const Fh = hM.rhs_form;
+
+       SECTION("dt_H")
+       {
+           auto const computed = At;
+           Eigen::MatrixXd expected(5,5);
+           expected <<
+                     -0.12575, 0.19111, -0.050545, -0.011227, -0.0035868,
+                     -0.70825, 0.99723, -0.1564, -0.10635, -0.026224,
+                     -0.41927, -0.1564, 0.89087, -0.20885, -0.10635,
+                     -0.2867, -0.10635, -0.20885, 0.78452, -0.18263,
+                     -0.13024, -0.026224, -0.10635, -0.18263, 0.44544;
+           auto const diff = computed - expected;
+           auto const diff_norm = diff.lpNorm<Eigen::Infinity>();
+              REQUIRE_THAT(diff_norm, Catch::Matchers::WithinAbs(0., 1e-4));
+       }SECTION("Hdt_dt")
+       {
+           auto const computed = Wt;
+           Eigen::MatrixXd expected(5,5);
+           expected <<
+                     0.89936, 0.36873, -0.62389, -0.11543, -0.52877,
+                     -1.43, 0.43738, 1.5011, -0.39304, -0.11543,
+                     0.43738, -1.3589, 0.04434, 1.5011, -0.62389,
+                     0.071088, 0.48172, -1.3589, 0.43738, 0.36873,
+                     0.02217, 0.071088, 0.43738, -1.43, 0.89936;
+           auto const diff = computed - expected;
+           auto const diff_norm = diff.lpNorm<Eigen::Infinity>();
+           REQUIRE_THAT(diff_norm, Catch::Matchers::WithinAbs(0., 1e-4));
+       }SECTION("I_H")
+       {
+           auto const computed = Mt;
+           Eigen::MatrixXd expected(5,5);
+           expected <<
+                     0.094731, -0.053269, -0.015677, -0.004206, -0.00058101,
+                     0.20271, 0.097697, -0.13247, -0.033489, -0.004206,
+                     0.097697, 0.27294, 0.064208, -0.13247, -0.015677,
+                     0.070234, 0.16191, 0.27294, 0.097697, -0.053269,
+                     0.032104, 0.070234, 0.097697, 0.20271, 0.094731;
+           auto const diff = computed - expected;
+           auto const diff_norm = diff.lpNorm<Eigen::Infinity>();
+           REQUIRE_THAT(diff_norm, Catch::Matchers::WithinAbs(0., 1e-4));
+       }SECTION("I0_H1(constant, tested with linear)")
+       {
+           auto const computed = Fh;
+           Eigen::MatrixXd expected(5,4);
+           expected <<
+                 0.062877, -0.032679, -0.0074071, -0.0017934,
+                 0.35412, -0.14449, -0.066288, -0.013112,
+                 0.20964, 0.28784, -0.1576, -0.053176,
+                 0.14335, 0.19652, 0.30095, -0.091313,
+                 0.065118, 0.07823, 0.13141, 0.22272;
+           auto const diff = computed - expected;
+           auto const diff_norm = diff.lpNorm<Eigen::Infinity>();
+           REQUIRE_THAT(diff_norm, Catch::Matchers::WithinAbs(0., 1e-4));
+       }
+    }
+    SECTION("Spline")
+    {
+        SECTION("dt_H")
+        {
+            auto const computed = modified_hilbert_transformation::get_dt_H_spline(tpoints);
+            Eigen::MatrixXd expected(6,6);
+            expected <<
+                     -0.21838, 0.28708, -0.053847, -0.010724, -0.0029812, -0.0011457,
+                     -0.6038, 0.61596, 0.11575, -0.098391, -0.021756, -0.0077663,
+                     -0.57947, 0.11575, 0.63332, -0.012163, -0.1233, -0.034138,
+                     -0.32364, -0.098391, -0.012163, 0.57427, 0.015599, -0.15568,
+                     -0.17947, -0.021756, -0.1233, 0.015599, 0.4042, -0.09528,
+                     -0.085134, -0.0077663, -0.034138, -0.15568, -0.09528, 0.37799;
+            auto const diff = computed - expected;
+            auto const diff_norm = diff.lpNorm<Eigen::Infinity>();
+            REQUIRE_THAT(diff_norm, Catch::Matchers::WithinAbs(0., 1e-4));
+        }SECTION("Hdt_dt")
+        {
+            auto const computed = modified_hilbert_transformation::get_Hdt_dt_spline(tpoints);
+            Eigen::MatrixXd expected(6,6);
+            expected <<
+                     1.5157, 0.10597, -0.84008, -0.21971, -0.048206, -0.51366,
+                     -1.9418, 0.71098, 1.5411, -0.22444, -0.037544, -0.048206,
+                     0.30074, -1.2214, 0.08575, 1.2791, -0.22444, -0.21971,
+                     0.091769, 0.30417, -1.1827, 0.08575, 1.5411, -0.84008,
+                     0.024352, 0.075956, 0.30417, -1.2214, 0.71098, 0.10597,
+                     0.0092962, 0.024352, 0.091769, 0.30074, -1.9418, 1.5157;
+            auto const diff = computed - expected;
+            auto const diff_norm = diff.lpNorm<Eigen::Infinity>();
+            REQUIRE_THAT(diff_norm, Catch::Matchers::WithinAbs(0., 1e-4));
+        }SECTION("I_H")
+        {
+            auto const computed = modified_hilbert_transformation::get_I_H_spline(tpoints);
+            Eigen::MatrixXd expected(6,6);
+            expected <<
+                     0.059545, -0.023465, -0.020384, -0.0041282, -0.00095459, -0.00014123,
+                     0.11635, 0.058125, -0.06865, -0.034004, -0.0066579, -0.00095459,
+                     0.094204, 0.1861, 0.071821, -0.10931, -0.034004, -0.0041282,
+                     0.053188, 0.12446, 0.23964, 0.071821, -0.06865, -0.020384,
+                     0.029709, 0.063454, 0.12446, 0.1861, 0.058125, -0.023465,
+                     0.014112, 0.029709, 0.053188, 0.094204, 0.11635, 0.059545;
+            auto const diff = computed - expected;
+            auto const diff_norm = diff.lpNorm<Eigen::Infinity>();
+            REQUIRE_THAT(diff_norm, Catch::Matchers::WithinAbs(0., 1e-4));
+        }SECTION("I1_H2 (p.w linear, tested with spline")
+        {
+            auto const computed = modified_hilbert_transformation::get_I_lin_H_spline(tpoints);
+            Eigen::MatrixXd expected(6,5);
+            expected <<
+                 0.054595, -0.034349, -0.0074254, -0.0020634, -0.00028642,
+                 0.15095, -0.0060813, -0.063957, -0.014761, -0.0019416,
+                 0.14487, 0.23186, -0.084799, -0.078718, -0.0085345,
+                 0.080911, 0.21102, 0.2171, -0.070038, -0.038919,
+                 0.044866, 0.10061, 0.16226, 0.15446, -0.02382,
+                 0.021283, 0.04645, 0.063519, 0.14136, 0.094498;
+            auto const diff = computed - expected;
+            auto const diff_norm = diff.lpNorm<Eigen::Infinity>();
+            REQUIRE_THAT(diff_norm, Catch::Matchers::WithinAbs(0., 1e-4));
+        }
+
+    }
+    REQUIRE(1==1);
+}
\ No newline at end of file

tests/z_matrices_tests.cpp

+//
+// Created by mreichelt on 23.05.25.
+//
+#include <catch2/catch_test_macros.hpp>
+#include <catch2/matchers/catch_matchers_floating_point.hpp>
+#include <catch2/benchmark/catch_benchmark_all.hpp>
+
+#include "modified_hilbert_transformation.hpp"
+#include "../submodules/Catch2/src/catch2/catch_test_macros.hpp"
+
+#include <iostream>
+
+template <typename T>
+std::vector<T> linspace(T a, T b, size_t N) {
+    T h = (b - a) / static_cast<T>(N-1);
+    std::vector<T> xs(N);
+    typename std::vector<T>::iterator x;
+    T val;
+    for (x = xs.begin(), val = a; x != xs.end(); ++x, val += h)
+        *x = val;
+    return xs;
+}
+
+TEST_CASE("Z-Matrices", "[hilbert_transformation][zMatrices]")
+{
+    //test against matlab produced matrices
+    auto const tpoints = linspace(0.,2., 5);
+    SECTION("Z0")
+    {
+        auto const Z0 = modified_hilbert_transformation::get_Z_0(tpoints);
+
+        Eigen::MatrixXd Z0_ref(5,4);
+        Z0_ref <<
+             1, 0, 0, 0,
+             0, 2, 0, 0,
+             0, -1, 3, 0,
+             0, 0, -2, 4,
+             0, 0, 0, -3;
+        auto const Diff = Z0 - Z0_ref;
+        auto const diff_norm = Diff.lpNorm<Eigen::Infinity>();
+        REQUIRE_THAT(diff_norm, Catch::Matchers::WithinAbs(0., 1e-10));
+    }
+
+    SECTION("Z1")
+    {
+        auto const Z1 = modified_hilbert_transformation::get_Z_1(tpoints);
+
+        Eigen::MatrixXd Z1_ref(5,4);
+        Z1_ref <<
+             -2, 0, 0, 0,
+             2, -2, 0, 0,
+             0, 2, -2, 0,
+             0, 0, 2, -2,
+             0, 0, 0, 2;
+        auto const Diff = Z1 - Z1_ref;
+        auto const diff_norm = Diff.lpNorm<Eigen::Infinity>();
+        REQUIRE_THAT(diff_norm, Catch::Matchers::WithinAbs(0., 1e-10));
+    }
+
+    SECTION("Z20")
+    {
+        auto const Z20 = modified_hilbert_transformation::get_Z_20(tpoints);
+
+        Eigen::MatrixXd Z20_ref(6,5);
+        Z20_ref <<
+                 1, 0, 0, 0, 0,
+                 0, 1, 0, 0, 0,
+                 0, 0, 1.5, 0, 0,
+                 0, 0, -0.5, 2, 0,
+                 0, 0, 0, -1, 4,
+                 0, 0, 0, 0, -3;
+        auto const Diff = Z20 - Z20_ref;
+        auto const diff_norm = Diff.lpNorm<Eigen::Infinity>();
+        REQUIRE_THAT(diff_norm, Catch::Matchers::WithinAbs(0., 1e-10));
+    }
+
+    SECTION("Z21")
+    {
+        auto const Z21 = modified_hilbert_transformation::get_Z_21(tpoints);
+
+        Eigen::MatrixXd Z21_ref(6,5);
+        Z21_ref <<
+             -2, 0, 0, 0, 0,
+             2, -1, 0, 0, 0,
+             0, 1, -1, 0, 0,
+             0, 0, 1, -1, 0,
+             0, 0, 0, 1, -2,
+             0, 0, 0, 0, 2;
+        auto const Diff = Z21 - Z21_ref;
+        auto const diff_norm = Diff.lpNorm<Eigen::Infinity>();
+        REQUIRE_THAT(diff_norm, Catch::Matchers::WithinAbs(0., 1e-10));
+    }
+
+    REQUIRE(1==1);
+}
\ No newline at end of file