Commit 94d352a4 authored by Timo Koch's avatar Timo Koch
Browse files

[common] Add monotone cubic spline and test

parent 0ab21de3
......@@ -24,6 +24,7 @@ intersectionmapper.hh
intrange.hh
loggingparametertree.hh
math.hh
monotonecubicspline.hh
numericdifferentiation.hh
optional.hh
parameters.hh
......
// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
// vi: set et ts=4 sw=4 sts=4:
/*****************************************************************************
* See the file COPYING for full copying permissions. *
* *
* This program is free software: you can redistribute it and/or modify *
* it under the terms of the GNU General Public License as published by *
* the Free Software Foundation, either version 3 of the License, or *
* (at your option) any later version. *
* *
* This program is distributed in the hope that it will be useful, *
* but WITHOUT ANY WARRANTY; without even the implied warranty of *
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
* GNU General Public License for more details. *
* *
* You should have received a copy of the GNU General Public License *
* along with this program. If not, see <http://www.gnu.org/licenses/>. *
*****************************************************************************/
/*!
* \file
* \ingroup Common
* \brief A monotone cubic spline
*/
#ifndef DUMUX_COMMON_MONOTONE_CUBIC_SPLINE_HH
#define DUMUX_COMMON_MONOTONE_CUBIC_SPLINE_HH
#include <cmath>
#include <cassert>
#include <iterator>
#include <vector>
#include <algorithm>
#include <dune/common/float_cmp.hh>
// Hermite basis functions
#include <dumux/common/cubicsplinehermitebasis.hh>
namespace Dumux {
/*!
* \ingroup Common
* \brief A monotone cubic spline
* \note Contruction 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>
class MonotoneCubicSpline
{
using Basis = CubicSplineHermiteBasis<Scalar>;
public:
/*!
* \brief Default constructor
*/
MonotoneCubicSpline() = default;
/*!
* \brief Contruct 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)
{
// check some requirements
assert (x.size() == y.size());
assert (x.size() >=2);
assert (std::is_sorted(x.begin(), x.end()));
updatePoints(x, y);
}
/*!
* \brief Create a monotone cubic spline from the control points (x[i], y[i])
* \param x a vector of x-coordinates
* \param y a vector of y-coordinates
*/
void updatePoints(const std::vector<Scalar>& x, const std::vector<Scalar>& y)
{
// save a copy of the control points
x_ = x;
y_ = y;
// the number of control points
numPoints_ = x.size();
// the slope at every control point
m_.resize(numPoints_);
// compute slopes (see Fritsch & Butland (1984), Eq. (5))
Scalar deltaX = (x[1]-x[0]);
Scalar secant = m_.front() = (y[1]-y[0])/deltaX;
Scalar prevDeltaX = deltaX;
Scalar prevSecant = secant;
for (int i = 1; i < numPoints_-1; ++i, prevSecant = secant, prevDeltaX = deltaX)
{
deltaX = (x[i+1]-x[i]);
secant = (y[i+1]-y[i])/deltaX;
const auto alpha = (prevDeltaX + 2*deltaX)/(3*(prevDeltaX + deltaX));
m_[i] = prevSecant*secant > 0.0 ? prevSecant*secant/(alpha*secant + (1.0-alpha)*prevSecant) : 0.0;
}
m_.back() = secant;
}
/*!
* \brief Evaluate the y value at a given x value
* \param x the x-coordinate
* \note We extrapolate linearly if out of bounds
*/
Scalar eval(const Scalar x) const
{
if (x <= x_.front())
return y_.front() + m_.front()*(x - x_.front());
else if (x > x_.back())
return y_.back() + m_.back()*(x - x_.back());
else
{
const auto lookUpIndex = std::distance(x_.begin(), std::lower_bound(x_.begin(), x_.end(), x));
assert(lookUpIndex != 0);
// 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_[lookUpIndex-1]*Basis::h00(t) + h*m_[lookUpIndex-1]*Basis::h10(t)
+ y_[lookUpIndex]*Basis::h01(t) + h*m_[lookUpIndex]*Basis::h11(t);
}
}
/*!
* \brief Evaluate the first derivative dy/dx at a given x value
* \param x the x-coordinate
* \note We extrapolate linearly if out of bounds
*/
Scalar evalDerivative(const Scalar x) const
{
if (x <= x_.front())
return m_.front();
else if (x > x_.back())
return m_.back();
else
{
const auto lookUpIndex = std::distance(x_.begin(), std::lower_bound(x_.begin(), x_.end(), x));
assert(lookUpIndex != 0);
// interpolate parametrization parameter t in [0,1]
const auto h = (x_[lookUpIndex] - x_[lookUpIndex-1]);
const auto t = (x - x_[lookUpIndex-1])/h;
const auto dtdx = 1.0/h;
return y_[lookUpIndex-1]*Basis::dh00(t)*dtdx + m_[lookUpIndex-1]*Basis::dh10(t)
+ y_[lookUpIndex]*Basis::dh01(t)*dtdx + m_[lookUpIndex]*Basis::dh11(t);
}
}
private:
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
};
} // end namespace Dumux
#endif
......@@ -2,3 +2,5 @@ dune_add_test(SOURCES test_spline.cc
LABELS unit)
dune_add_test(SOURCES test_cubicspline.cc
LABELS unit)
dune_add_test(SOURCES test_monotonecubicspline.cc
LABELS unit)
......@@ -23,6 +23,7 @@
#include <config.h>
#include <cmath>
#include <vector>
#include <numeric>
#include <algorithm>
#include <functional>
......
// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
// vi: set et ts=4 sw=4 sts=4:
/*****************************************************************************
* See the file COPYING for full copying permissions. *
* *
* This program is free software: you can redistribute it and/or modify *
* it under the terms of the GNU General Public License as published by *
* the Free Software Foundation, either version 3 of the License, or *
* (at your option) any later version. *
* *
* This program is distributed in the hope that it will be useful, *
* but WITHOUT ANY WARRANTY; without even the implied warranty of *
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
* GNU General Public License for more details. *
* *
* You should have received a copy of the GNU General Public License *
* along with this program. If not, see <http://www.gnu.org/licenses/>. *
*****************************************************************************/
/*!
* \file
* \brief Test the simple cubic spline implementation
*/
#include <config.h>
#include <cmath>
#include <vector>
#include <numeric>
#include <algorithm>
#include <functional>
#include <dune/common/exceptions.hh>
#include <dune/common/parallel/mpihelper.hh>
#include <dumux/common/monotonecubicspline.hh>
#include <dumux/io/gnuplotinterface.hh>
std::vector<double> linspace(const double begin, const double end, const double samples)
{
const double delta = (end-begin)/static_cast<double>(samples-1);
std::vector<double> vec(samples);
for (int i = 0; i < samples; ++i)
vec[i] = begin + i*delta;
return vec;
}
template<class Function>
std::vector<double> eval(const Function& f, const std::vector<double>& x)
{
auto y = x;
std::transform(x.begin(), x.end(), y.begin(), [&](const double x) { return f(x); });
return y;
}
int main(int argc, char** argv)
{
Dune::MPIHelper::instance(argc, argv);
// we test the spline interpolation against a sample function
const auto f = [](double x){ return x*x*x; };
const auto df = [](double x){ return 3*x*x; };
// create some test samples
const auto testPoints = linspace(0.0, 4.0, 1000);
const auto ref = eval(f, testPoints);
const auto refDeriv = eval(df, testPoints);
// create the spline sample points
const auto samplePoints = linspace(0.0, 5.0, 10);
const auto y = eval(f, samplePoints);
// create the spline
Dumux::MonotoneCubicSpline spline(samplePoints, y);
// evaluate spline and derivative
const auto result = eval([&](const double x) { return spline.eval(x); }, testPoints);
const auto resultDeriv = eval([&](const double x) { return spline.evalDerivative(x); }, testPoints);
// compute largest difference
auto diff = result; auto diffDeriv = result;
std::transform(result.begin(), result.end(), ref.begin(), diff.begin(), [](auto a, auto b){ return std::abs(a-b); });
std::transform(resultDeriv.begin(), resultDeriv.end(), refDeriv.begin(), diffDeriv.begin(), [](auto a, auto b){ return std::abs(a-b); });
const auto maxNorm = std::accumulate(diff.begin(), diff.end(), diff[0], [](auto a, auto b){ return std::max(a, b); })
/(*std::max_element(ref.begin(), ref.end()));
const auto maxNormDeriv = std::accumulate(diffDeriv.begin(), diffDeriv.end(), diffDeriv[0], [](auto a, auto b){ return std::max(a, b); })
/(*std::max_element(refDeriv.begin(), refDeriv.end()));
std::cout << "Maximum error: " << maxNorm << "\n";
std::cout << "Maximum error in derivative: " << maxNormDeriv << "\n";
if (maxNorm > 0.0008 || maxNormDeriv > 0.013)
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 = linspace(-1.0, 5.0, 1000);
const auto refPlot = eval(f, plotPoints);
const auto refDerivPlot = eval(df, plotPoints);
const auto resultPlot = eval([&](const double x) { return spline.eval(x); }, plotPoints);
const auto resultDerivPlot = eval([&](const double x) { return spline.evalDerivative(x); }, plotPoints);
Dumux::GnuplotInterface<double> gnuplot(false);
gnuplot.addDataSetToPlot(plotPoints, refPlot, "exp_reference");
gnuplot.addDataSetToPlot(plotPoints, refDerivPlot, "exp_reference_derivative");
gnuplot.addDataSetToPlot(plotPoints, resultPlot, "monotspline");
gnuplot.addDataSetToPlot(plotPoints, resultDerivPlot, "monotspline_derivative");
gnuplot.plot("monotspline");
return 0;
}
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment