diff --git a/tests/test_FunctionEvaluation.cpp b/tests/test_FunctionEvaluation.cpp index 88c9e552af82fdfecc1d2da1a78cbefd2ffae969..7893ba07eec45c4b33c842454172dde6446e46ac 100644 --- a/tests/test_FunctionEvaluation.cpp +++ b/tests/test_FunctionEvaluation.cpp @@ -38,16 +38,29 @@ TEST_CASE("QuadratureFormula", "[analysis]") auto point_list = quadrature_formula.pointList(); auto weight_list = quadrature_formula.weightList(); - auto x = [&a, &b](auto x_hat) { return 0.5 * (b - a) * (x_hat + 1); }; + double alpha = 0.5 * (b - a); + double beta = 0.5 * (a + b); + + auto x = [&alpha, &beta](auto x_hat) { return alpha * x_hat + beta; }; auto value = weight_list[0] * f(x(point_list[0])); for (size_t i = 1; i < weight_list.size(); ++i) { value += weight_list[i] * f(x(point_list[i])); } - return 0.5 * (b - a) * value; + + return alpha * value; + }; + + auto get_order = [&integrate](auto f, auto quadrature_formula, const double a, const double b, + const double exact_value) { + double int_ab = integrate(f, quadrature_formula, a, b); + double int_first_half = integrate(f, quadrature_formula, a, 0.5 * (a + b)); + double int_second_half = integrate(f, quadrature_formula, 0.5 * (a + b), b); + + return -std::log((int_first_half + int_second_half - exact_value) / (int_ab - exact_value)) / std::log(2); }; - auto p0 = [](const TinyVector<1>&) { return 3; }; + auto p0 = [](const TinyVector<1>&) { return 1; }; auto p1 = [](const TinyVector<1>& X) { const double x = X[0]; return 2 * x + 1; @@ -56,6 +69,14 @@ TEST_CASE("QuadratureFormula", "[analysis]") const double x = X[0]; return 3 * x * x + 2 * x + 1; }; + auto p3 = [](const TinyVector<1>& X) { + const double x = X[0]; + return 4 * std::pow(x, 3) + 3 * x * x + 2 * x + 1; + }; + auto p4 = [](const TinyVector<1>& X) { + const double x = X[0]; + return 5 * std::pow(x, 4) + 4 * std::pow(x, 3) + 3 * x * x + 2 * x + 1; + }; SECTION("exact integration") { @@ -67,8 +88,11 @@ TEST_CASE("QuadratureFormula", "[analysis]") REQUIRE(l0.numberOfPoints() == 1); REQUIRE(l1.numberOfPoints() == 1); - REQUIRE(integrate(p0, l1, 0.5, 1) == Catch::Approx(1.5)); + REQUIRE(integrate(p0, l1, 0.5, 1) == Catch::Approx(0.5)); REQUIRE(integrate(p1, l1, 0, 1) == Catch::Approx(2)); + REQUIRE(integrate(p2, l1, 0, 1) != Catch::Approx(3)); + + REQUIRE(get_order(p2, l1, 0, 1, 3) == Catch::Approx(2)); } SECTION("order 2 and 3") @@ -79,9 +103,13 @@ TEST_CASE("QuadratureFormula", "[analysis]") REQUIRE(l2.numberOfPoints() == 2); REQUIRE(l3.numberOfPoints() == 2); - REQUIRE(integrate(p0, l3, 0.5, 1) == Catch::Approx(1.5)); + REQUIRE(integrate(p0, l3, 0.5, 1) == Catch::Approx(0.5)); REQUIRE(integrate(p1, l3, 0, 1) == Catch::Approx(2)); REQUIRE(integrate(p2, l3, 0, 1) == Catch::Approx(3)); + REQUIRE(integrate(p3, l3, 0, 1) == Catch::Approx(4)); + REQUIRE(integrate(p4, l3, 0, 1) != Catch::Approx(5)); + + REQUIRE(get_order(p4, l3, 0, 1, 5) == Catch::Approx(4)); } } }