Commit 2df2f565 by Kilian Weishaupt

### Merge branch 'feature/monotone-spline-inverse' into 'master'

```[common] Add inverse of monotone cubic spline

See merge request !2250```
parents 0626ef73 a23766d9
 ... ... @@ -35,12 +35,15 @@ // Hermite basis functions #include // for inversion #include 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 ... ... @@ -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& x, const std::vector y) MonotoneCubicSpline(const std::vector& x, const std::vector& 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& x, const std::vector& 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 x_; //!< the x-coordinates std::vector y_; //!< the y-coordinates std::vector 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 ... ...
 ... ... @@ -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); ... ...
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!