Merge pull request #18 from i05nagai/add-extrapolation

Add extrapolation

recipe/quadrature/richardson_extrapolation.cc

@@ -0,0 +1,46 @@
+#include "recipe/quadrature/richardson_extrapolation.h"
+#include <cassert>
+#include <cmath>
+#include "recipe/quadrature/quadrature.h"
+
+namespace recipe {
+namespace quadrature {
+
+double RichardsonExtrapolation(QuadratureRule quadrature_rule,
+                               double (*integrand)(double),
+                               double left_endpoint, double right_endpoint,
+                               unsigned int num_of_partition) {
+  assert(left_endpoint < right_endpoint);
+
+  double approximator_of_half_step_size =
+      0;  // richardson extrapolation uses the approximation value of a half
+          // step size.
+  double approximator_of_normal_step_size =
+      0;  // the approximation value of a normal step size.
+  unsigned int approx_order = 0;
+
+  switch (quadrature_rule) {
+    case Trapezoidal: {
+      approx_order = 2;
+      approximator_of_normal_step_size = recipe::quadrature::TrapezoidalRule(
+          integrand, left_endpoint, right_endpoint, num_of_partition);
+      approximator_of_half_step_size = recipe::quadrature::TrapezoidalRule(
+          integrand, left_endpoint, right_endpoint, 2 * num_of_partition);
+      break;
+    }
+    case Simpson: {
+      approx_order = 4;
+      approximator_of_normal_step_size = recipe::quadrature::SimpsonRule(
+          integrand, left_endpoint, right_endpoint, num_of_partition);
+      approximator_of_half_step_size = recipe::quadrature::SimpsonRule(
+          integrand, left_endpoint, right_endpoint, 2 * num_of_partition);
+      break;
+    }
+    default: { assert("An undefined quadrature rule was specified."); }
+  }
+  return (std::pow(2.0, approx_order) * approximator_of_half_step_size -
+          approximator_of_normal_step_size) /
+         (pow(2.0, approx_order) - 1.0);
+};
+}  // namespace quadrature
+}  // namespace recipe

recipe/quadrature/richardson_extrapolation.h

@@ -0,0 +1,17 @@
+namespace recipe {
+namespace quadrature {
+/*! \enum QuadratureRule
+ *
+ *  QuadratureRule enumerator lists all the quadrature rules that have been
+ * implemented. Avairable are: Trapezoidal, Simpson. Please put a name in the
+ * definition when you implements a new quadrature rule.
+ *
+ */
+enum QuadratureRule { Trapezoidal, Simpson };
+
+double RichardsonExtrapolation(QuadratureRule quadrature_rule,
+                               double (*integrand)(double),
+                               double left_endpoint, double right_endpoint,
+                               unsigned int num_of_partition);
+}  // namespace quadrature
+}  // namespace recipe

recipe/quadrature/richardson_extrapolation_test.cc

@@ -0,0 +1,82 @@
+#include "recipe/quadrature/richardson_extrapolation.h"
+#include <gtest/gtest.h>
+#include "recipe/test_util/test_util.h"
+
+namespace recipe {
+namespace quadrature {
+
+double MonomOf6thDegree(double argument) {
+  return argument * argument * argument * argument * argument * argument;
+}  // f(x)=x^6
+double MonomOf8thDegree(double argument) {
+  return argument * argument * argument * argument * argument * argument *
+         argument * argument;
+}  // f(x)=x^8
+
+#ifndef NDEBUG
+TEST(RichardsonExtrapolation, AssertIntervalForSimpson) {
+  double left_endpoint = 1.0;
+  double right_endpoint = 0;
+  unsigned int num_of_partition = 100;
+
+  EXPECT_ASSERT_FAIL(recipe::quadrature::RichardsonExtrapolation(
+      recipe::quadrature::QuadratureRule::Simpson, &MonomOf6thDegree,
+      left_endpoint, right_endpoint, num_of_partition));
+}
+#endif
+
+#ifndef NDEBUG
+TEST(RichardsonExtrapolation, AssertIntervalForTrapezoidal) {
+  double left_endpoint = 1.0;
+  double right_endpoint = 0;
+  unsigned int num_of_partition = 100;
+
+  EXPECT_ASSERT_FAIL(recipe::quadrature::RichardsonExtrapolation(
+      recipe::quadrature::QuadratureRule::Trapezoidal, &MonomOf6thDegree,
+      left_endpoint, right_endpoint, num_of_partition));
+}
+#endif
+
+// test for extrapolation with trapezoidal rule on the integrations of the
+// monomial of 6th degrees from 0 to 1.
+TEST(RichardsonExtrTest, TrapezoidalWithMonomOf6thDegreeAsIntegrand) {
+  unsigned int num_of_partition =
+      100;                      // denoted by N in the following comments.
+  double left_endpoint = 0;     // integraton from zero makes it easy to compute
+                                // the closed form for polynomial integrands.
+  double right_endpoint = 1.0;  // denoted by a.
+
+  double approximator = recipe::quadrature::RichardsonExtrapolation(
+      recipe::quadrature::QuadratureRule::Trapezoidal, &MonomOf6thDegree,
+      left_endpoint, right_endpoint, num_of_partition);
+
+  // error turns out to be bounded by a^7/(24N^4).
+  double error_boundary = 1.0 / (24.0 * num_of_partition * num_of_partition *
+                                 num_of_partition * num_of_partition);
+
+  EXPECT_NEAR(approximator, 1.0 / 7.0, error_boundary);
+}
+
+// test for extrapolation with Simpson rule on the integrations of the monomial
+// of 8th degrees from 0 to 1.
+TEST(RichardsonExtrTest, SimpsonWithMonomOf8thDegreeAsIntegrand) {
+  unsigned int num_of_partition =
+      100;                      // denoted by N in the following comments.
+  double left_endpoint = 0;     // integraton from zero makes it easy to compute
+                                // the closed form for polynomial integrands.
+  double right_endpoint = 1.0;  // denoted by a.
+
+  double approximator = recipe::quadrature::RichardsonExtrapolation(
+      recipe::quadrature::QuadratureRule::Simpson, &MonomOf8thDegree,
+      left_endpoint, right_endpoint, num_of_partition);
+
+  // error turns out to be bounded by 5a^9/(95N^6).
+  double error_boundary =
+      5.0 / (95.0 * num_of_partition * num_of_partition * num_of_partition *
+             num_of_partition * num_of_partition * num_of_partition);
+
+  EXPECT_NEAR(approximator, 1.0 / 9.0, error_boundary);
+}
+
+}  // namespace quadrature
+}  // namespace recipe