diff --git a/dumux/common/monotonecubicspline.hh b/dumux/common/monotonecubicspline.hh
index 8a6bcdc7a32edd02303e780a6d3cc270bae7fbcd..46a8727feae7921501989c7283dc09a47af6f781 100644
--- a/dumux/common/monotonecubicspline.hh
+++ b/dumux/common/monotonecubicspline.hh
@@ -35,12 +35,15 @@
// Hermite basis functions
#include <dumux/common/cubicsplinehermitebasis.hh>
+// for inversion
+#include <dumux/nonlinear/findscalarroot.hh>
+
namespace Dumux {
/*!
* \ingroup Common
* \brief A monotone cubic spline
- * \note Contruction after Fritsch & Butland (1984) (see https://doi.org/10.1137/0905021)
+ * \note Construction after Fritsch & Butland (1984) (see https://doi.org/10.1137/0905021)
* \note The resulting interpolation is globally monotone but only C^1
*/
template<class Scalar = double>
@@ -54,18 +57,13 @@ public:
MonotoneCubicSpline() = default;
/*!
- * \brief Contruct a monotone cubic spline from the control points (x[i], y[i])
+ * \brief Construct a monotone cubic spline from the control points (x[i], y[i])
* \note if the data set is monotone, monotonicity is preserved
* \param x a vector of x-coordinates
* \param y a vector of y-coordinates
*/
- MonotoneCubicSpline(const std::vector<Scalar>& x, const std::vector<Scalar> y)
+ MonotoneCubicSpline(const std::vector<Scalar>& x, const std::vector<Scalar>& y)
{
- // check some requirements
- assert (x.size() == y.size());
- assert (x.size() >=2);
- assert (std::is_sorted(x.begin(), x.end()));
-
updatePoints(x, y);
}
@@ -76,6 +74,11 @@ public:
*/
void updatePoints(const std::vector<Scalar>& x, const std::vector<Scalar>& y)
{
+ // check some requirements
+ assert (x.size() == y.size());
+ assert (x.size() >=2);
+ assert (std::is_sorted(x.begin(), x.end()));
+
// save a copy of the control points
x_ = x;
y_ = y;
@@ -83,6 +86,9 @@ public:
// the number of control points
numPoints_ = x.size();
+ // whether we are increasing
+ increasing_ = y_.back() > y_.front();
+
// the slope at every control point
m_.resize(numPoints_);
@@ -150,11 +156,66 @@ public:
}
}
+ /*!
+ * \brief Evaluate the inverse function
+ * \param x the x-coordinate
+ * \note We extrapolate linearly if out of bounds
+ * \note Throws exception if inverse could not be found (e.g. not unique)
+ */
+ Scalar evalInverse(const Scalar y) const
+ {
+ if (increasing_)
+ {
+ if (y <= y_.front())
+ return x_.front() + (y - y_.front())/m_.front();
+ else if (y > y_.back())
+ return x_.back() + (y - y_.back())/m_.back();
+ else
+ {
+ const auto lookUpIndex = std::distance(y_.begin(), std::lower_bound(y_.begin(), y_.end(), y));
+ assert(lookUpIndex != 0);
+
+ return evalInverse_(y, lookUpIndex);
+ }
+ }
+
+ else
+ {
+ if (y >= y_.front())
+ return x_.front() + (y - y_.front())/m_.front();
+ else if (y < y_.back())
+ return x_.back() + (y - y_.back())/m_.back();
+ else
+ {
+ const auto lookUpIndex = y_.size() - std::distance(y_.rbegin(), std::lower_bound(y_.rbegin(), y_.rend(), y));
+ assert(lookUpIndex != 0);
+
+ return evalInverse_(y, lookUpIndex);
+ }
+ }
+ }
+
private:
+ Scalar evalInverse_(const Scalar y, const std::size_t lookUpIndex) const
+ {
+ auto localPolynomial = [&](const auto x) {
+ // interpolate parametrization parameter t in [0,1]
+ const auto h = (x_[lookUpIndex] - x_[lookUpIndex-1]);
+ const auto t = (x - x_[lookUpIndex-1])/h;
+ return y - (y_[lookUpIndex-1]*Basis::h00(t) + h*m_[lookUpIndex-1]*Basis::h10(t)
+ + y_[lookUpIndex]*Basis::h01(t) + h*m_[lookUpIndex]*Basis::h11(t));
+ };
+
+ // use an epsilon for the bracket
+ const auto eps = (x_[lookUpIndex]-x_[lookUpIndex-1])*1e-5;
+ return findScalarRootBrent(x_[lookUpIndex-1]-eps, x_[lookUpIndex]+eps, localPolynomial);
+ }
+
std::vector<Scalar> x_; //!< the x-coordinates
std::vector<Scalar> y_; //!< the y-coordinates
std::vector<Scalar> m_; //!< the slope for each control point
std::size_t numPoints_; //!< the number of control points
+ bool increasing_; //!< if we are increasing monotone or not
};
} // end namespace Dumux
diff --git a/test/common/spline/test_monotonecubicspline.cc b/test/common/spline/test_monotonecubicspline.cc
index 378efe26487699a2b0c10a81d6b2616e8ea87b82..4ff362388145eb75db28ae7a477208a619ff54b7 100644
--- a/test/common/spline/test_monotonecubicspline.cc
+++ b/test/common/spline/test_monotonecubicspline.cc
@@ -79,6 +79,34 @@ int main(int argc, char** argv)
if (maxNorm > 0.0008 || maxNormDeriv > 0.013)
DUNE_THROW(Dune::Exception, "Maximum error in spline interpolation too large!");
+ // test inverse by evaluating (x = f^-1(f(x))) for monotonically increasing function
+ {
+ const auto resultX = eval([&](double x){ return spline.evalInverse(spline.eval(x)); }, testPoints);
+ auto diffInverse = resultX;
+ std::transform(resultX.begin(), resultX.end(), testPoints.begin(), diffInverse.begin(), [](auto a, auto b){ return std::abs(a-b); });
+ const auto maxNormInverse = std::accumulate(diffInverse.begin(), diffInverse.end(), diffInverse[0], [](auto a, auto b){ return std::max(a, b); })
+ /(*std::max_element(testPoints.begin(), testPoints.end()));
+
+
+ std::cout << "Maximum error in identity using the inverse (mon. incr.): " << std::scientific << maxNormInverse << "\n";
+ if (maxNormInverse > 1e-13)
+ DUNE_THROW(Dune::Exception, "Maximum error in spline interpolation too large!");
+ }
+ // test inverse by evaluating (x = f^-1(f(x))) for monotonically decreasing function
+ {
+ auto reverseTest = testPoints;
+ std::reverse(reverseTest.begin(), reverseTest.end());
+ const auto resultX = eval([&](double x){ return spline.evalInverse(spline.eval(x)); }, reverseTest);
+ auto diffInverse = resultX;
+ std::transform(resultX.begin(), resultX.end(), reverseTest.begin(), diffInverse.begin(), [](auto a, auto b){ return std::abs(a-b); });
+ const auto maxNormInverse = std::accumulate(diffInverse.begin(), diffInverse.end(), diffInverse[0], [](auto a, auto b){ return std::max(a, b); })
+ /(*std::max_element(reverseTest.begin(), reverseTest.end()));
+
+ std::cout << "Maximum error in identity using the inverse (mon. decr.): " << std::scientific << maxNormInverse << "\n";
+ if (maxNormInverse > 1e-13)
+ DUNE_THROW(Dune::Exception, "Maximum error in spline interpolation too large!");
+ }
+
// plot with Gnuplot (plot a bit more so we can see the linear extension)
const auto plotPoints = Dumux::linspace(-1.0, 5.0, 1000);
const auto refPlot = eval(f, plotPoints);