MonotoneCubicInterpolation Class Reference

This class interpolates values given a set of data pairs and an abscissa. More...

#include <MonotoneCubicInterpolation.h>

Public Member Functions
	MonotoneCubicInterpolation ()
	Empty constructor.
	MonotoneCubicInterpolation (const std::vector< Real > &x, const std::vector< Real > &y)
	Constructor, Takes two vectors of points for which to apply the fit.
virtual	~MonotoneCubicInterpolation ()=default
virtual void	setData (const std::vector< Real > &x, const std::vector< Real > &y)
	Method generally used when MonotoneCubicInterpolation object was created using the empty constructor.
virtual Real	sample (const Real &x) const
	This function will take an independent variable input and will return the dependent variable based on the generated fit.
virtual Real	sampleDerivative (const Real &x) const
	This function will take an independent variable input and will return the derivative of the dependent variable with respect to the independent variable based on the generated fit.
virtual Real	sample2ndDerivative (const Real &x) const
	This function will take an independent variable input and will return the second derivative of the dependent variable with respect to the independent variable based on the generated fit.
virtual void	dumpCSV (std::string filename, const std::vector< Real > &xnew)
	This function takes an array of independent variable values and writes a CSV file with values corresponding to y, y', and y''.
virtual unsigned int	getSampleSize ()
	This method returns the length of the independent variable vector.

Protected Member Functions
virtual void	errorCheck ()
Real	sign (const Real &x) const
Real	phi (const Real &t) const
Real	psi (const Real &t) const
Real	phiPrime (const Real &t) const
Real	psiPrime (const Real &t) const
Real	phiDoublePrime (const Real &t) const
Real	psiDoublePrime (const Real &t) const
Real	h1 (const Real &xhi, const Real &xlo, const Real &x) const
Real	h2 (const Real &xhi, const Real &xlo, const Real &x) const
Real	h3 (const Real &xhi, const Real &xlo, const Real &x) const
Real	h4 (const Real &xhi, const Real &xlo, const Real &x) const
Real	h1Prime (const Real &xhi, const Real &xlo, const Real &x) const
Real	h2Prime (const Real &xhi, const Real &xlo, const Real &x) const
Real	h3Prime (const Real &xhi, const Real &xlo, const Real &x) const
Real	h4Prime (const Real &xhi, const Real &xlo, const Real &x) const
Real	h1DoublePrime (const Real &xhi, const Real &xlo, const Real &x) const
Real	h2DoublePrime (const Real &xhi, const Real &xlo, const Real &x) const
Real	h3DoublePrime (const Real &xhi, const Real &xlo, const Real &x) const
Real	h4DoublePrime (const Real &xhi, const Real &xlo, const Real &x) const
virtual Real	p (const Real &xhi, const Real &xlo, const Real &fhi, const Real &flo, const Real &dhi, const Real &dlo, const Real &x) const
virtual Real	pPrime (const Real &xhi, const Real &xlo, const Real &fhi, const Real &flo, const Real &dhi, const Real &dlo, const Real &x) const
virtual Real	pDoublePrime (const Real &xhi, const Real &xlo, const Real &fhi, const Real &flo, const Real &dhi, const Real &dlo, const Real &x) const
virtual void	initialize_derivs ()
virtual void	modify_derivs (const Real &alpha, const Real &beta, const Real &delta, Real &yp_lo, Real &yp_hi)
virtual void	solve ()
virtual void	findInterval (const Real &x, unsigned int &klo, unsigned int &khi) const

Protected Attributes
std::vector< Real >	_x
std::vector< Real >	_y
std::vector< Real >	_h
std::vector< Real >	_yp
std::vector< Real >	_delta
std::vector< Real >	_alpha
std::vector< Real >	_beta
unsigned int	_n_knots
unsigned int	_n_intervals
unsigned int	_internal_knots

Detailed Description

This class interpolates values given a set of data pairs and an abscissa.

The interpolation is cubic with at least C1 continuity; C2 continuity can be violated in favor of ensuring the interpolation is monotonic. The algorithm used is laid out in Fritsch and Carlson, SIAM J. Numer. Anal. Vol. 17(2) April 1980

Definition at line 25 of file MonotoneCubicInterpolation.h.

Constructor & Destructor Documentation

MonotoneCubicInterpolation() [1/2]

MonotoneCubicInterpolation::MonotoneCubicInterpolation

(

)

Empty constructor.

Definition at line 20 of file MonotoneCubicInterpolation.C.

{}

MonotoneCubicInterpolation() [2/2]

MonotoneCubicInterpolation::MonotoneCubicInterpolation	(	const std::vector< Real > &	x,
		const std::vector< Real > &	y
	)

Constructor, Takes two vectors of points for which to apply the fit.

One should be of the independent variable while the other should be of the dependent variable. These values should have a one-to-one correspondence, e.g. the vectors must be of the same size.

Definition at line 22 of file MonotoneCubicInterpolation.C.

  : _x(x), _y(y)
{
  errorCheck();
  solve();
}

~MonotoneCubicInterpolation()

virtual MonotoneCubicInterpolation::~MonotoneCubicInterpolation ( )

virtualdefault

Member Function Documentation

dumpCSV()

void MonotoneCubicInterpolation::dumpCSV ( std::string  filename, const std::vector< Real > &  xnew  )

virtual

This function takes an array of independent variable values and writes a CSV file with values corresponding to y, y', and y''.

This can be used for sanity checks of the interpolation curve.

Definition at line 375 of file MonotoneCubicInterpolation.C.

{
  unsigned int n = xnew.size();
  std::vector<Real> ynew(n), ypnew(n), yppnew(n);
 
  std::ofstream out(filename.c_str());
  for (unsigned int i = 0; i < n; ++i)
  {
    ynew[i] = sample(xnew[i]);
    ypnew[i] = sampleDerivative(xnew[i]);
    yppnew[i] = sample2ndDerivative(xnew[i]);
    out << xnew[i] << ", " << ynew[i] << ", " << ypnew[i] << ", " << yppnew[i] << "\n";
  }
 
  out << std::flush;
  out.close();
}

errorCheck()

void MonotoneCubicInterpolation::errorCheck ( )

protectedvirtual

Definition at line 40 of file MonotoneCubicInterpolation.C.

{
  if (_x.size() != _y.size())
    throw std::domain_error("MonotoneCubicInterpolation: x and y vectors are not the same length");
 
  if (_x.size() < 3)
    throw std::domain_error("MonotoneCubicInterpolation: " + Moose::stringify(_x.size()) +
                            " points is not enough data for a cubic interpolation");
 
  bool error = false;
  for (unsigned i = 0; !error && i + 1 < _x.size(); ++i)
    if (_x[i] >= _x[i + 1])
      error = true;
 
  if (error)
    throw std::domain_error("x-values are not strictly increasing");
}

Referenced by MonotoneCubicInterpolation(), and setData().

findInterval()

void MonotoneCubicInterpolation::findInterval ( const Real &  x, unsigned int &  klo, unsigned int &  khi  ) const

protectedvirtual

Definition at line 330 of file MonotoneCubicInterpolation.C.

{
  klo = 0;
  khi = _n_knots - 1;
  while (khi - klo > 1)
  {
    unsigned int k = (khi + klo) >> 1;
    if (_x[k] > x)
      khi = k;
    else
      klo = k;
  }
}

Referenced by sample(), sample2ndDerivative(), and sampleDerivative().

getSampleSize()

unsigned int MonotoneCubicInterpolation::getSampleSize ( )

virtual

This method returns the length of the independent variable vector.

Definition at line 394 of file MonotoneCubicInterpolation.C.

{
  return _x.size();
}

h1()

Real MonotoneCubicInterpolation::h1 ( const Real &  xhi, const Real &  xlo, const Real &  x  ) const

protected

Definition at line 106 of file MonotoneCubicInterpolation.C.

{
  Real h = xhi - xlo;
  Real t = (xhi - x) / h;
  return phi(t);
}

Referenced by p().

h1DoublePrime()

Real MonotoneCubicInterpolation::h1DoublePrime ( const Real &  xhi, const Real &  xlo, const Real &  x  ) const

protected

Definition at line 123 of file MonotoneCubicInterpolation.C.

{
  Real h = xhi - xlo;
  Real t = (xhi - x) / h;
  Real tPrime = -1. / h;
  return phiDoublePrime(t) * tPrime * tPrime;
}

Referenced by pDoublePrime().

h1Prime()

Real MonotoneCubicInterpolation::h1Prime ( const Real &  xhi, const Real &  xlo, const Real &  x  ) const

protected

Definition at line 114 of file MonotoneCubicInterpolation.C.

{
  Real h = xhi - xlo;
  Real t = (xhi - x) / h;
  Real tPrime = -1. / h;
  return phiPrime(t) * tPrime;
}

Referenced by pPrime().

h2()

Real MonotoneCubicInterpolation::h2 ( const Real &  xhi, const Real &  xlo, const Real &  x  ) const

protected

Definition at line 132 of file MonotoneCubicInterpolation.C.

{
  Real h = xhi - xlo;
  Real t = (x - xlo) / h;
  return phi(t);
}

Referenced by p().

h2DoublePrime()

Real MonotoneCubicInterpolation::h2DoublePrime ( const Real &  xhi, const Real &  xlo, const Real &  x  ) const

protected

Definition at line 149 of file MonotoneCubicInterpolation.C.

{
  Real h = xhi - xlo;
  Real t = (x - xlo) / h;
  Real tPrime = 1. / h;
  return phiDoublePrime(t) * tPrime * tPrime;
}

Referenced by pDoublePrime().

h2Prime()

Real MonotoneCubicInterpolation::h2Prime ( const Real &  xhi, const Real &  xlo, const Real &  x  ) const

protected

Definition at line 140 of file MonotoneCubicInterpolation.C.

{
  Real h = xhi - xlo;
  Real t = (x - xlo) / h;
  Real tPrime = 1. / h;
  return phiPrime(t) * tPrime;
}

Referenced by pPrime().

h3()

Real MonotoneCubicInterpolation::h3 ( const Real &  xhi, const Real &  xlo, const Real &  x  ) const

protected

Definition at line 158 of file MonotoneCubicInterpolation.C.

{
  Real h = xhi - xlo;
  Real t = (xhi - x) / h;
  return -h * psi(t);
}

Referenced by p().

h3DoublePrime()

Real MonotoneCubicInterpolation::h3DoublePrime ( const Real &  xhi, const Real &  xlo, const Real &  x  ) const

protected

Definition at line 175 of file MonotoneCubicInterpolation.C.

{
  Real h = xhi - xlo;
  Real t = (xhi - x) / h;
  Real tPrime = -1. / h;
  return psiDoublePrime(t) * tPrime;
}

Referenced by pDoublePrime().

h3Prime()

Real MonotoneCubicInterpolation::h3Prime ( const Real &  xhi, const Real &  xlo, const Real &  x  ) const

protected

Definition at line 166 of file MonotoneCubicInterpolation.C.

{
  Real h = xhi - xlo;
  Real t = (xhi - x) / h;
  Real tPrime = -1. / h;
  return -h * psiPrime(t) * tPrime; // psiPrime(t)
}

Referenced by pPrime().

h4()

Real MonotoneCubicInterpolation::h4 ( const Real &  xhi, const Real &  xlo, const Real &  x  ) const

protected

Definition at line 184 of file MonotoneCubicInterpolation.C.

{
  Real h = xhi - xlo;
  Real t = (x - xlo) / h;
  return h * psi(t);
}

Referenced by p().

h4DoublePrime()

Real MonotoneCubicInterpolation::h4DoublePrime ( const Real &  xhi, const Real &  xlo, const Real &  x  ) const

protected

Definition at line 201 of file MonotoneCubicInterpolation.C.

{
  Real h = xhi - xlo;
  Real t = (x - xlo) / h;
  Real tPrime = 1. / h;
  return psiDoublePrime(t) * tPrime;
}

Referenced by pDoublePrime().

h4Prime()

Real MonotoneCubicInterpolation::h4Prime ( const Real &  xhi, const Real &  xlo, const Real &  x  ) const

protected

Definition at line 192 of file MonotoneCubicInterpolation.C.

{
  Real h = xhi - xlo;
  Real t = (x - xlo) / h;
  Real tPrime = 1. / h;
  return h * psiPrime(t) * tPrime; // psiPrime(t)
}

Referenced by pPrime().

initialize_derivs()

void MonotoneCubicInterpolation::initialize_derivs ( )

protectedvirtual

Definition at line 249 of file MonotoneCubicInterpolation.C.

{
  for (unsigned int i = 1; i < _n_knots - 1; ++i)
    _yp[i] = (std::pow(_h[i - 1], 2) * _y[i + 1] - std::pow(_h[i], 2) * _y[i - 1] -
              _y[i] * (_h[i - 1] - _h[i]) * (_h[i - 1] + _h[i])) /
             (_h[i - 1] * _h[i] * (_h[i - 1] * _h[i]));
 
  _yp[0] = (-std::pow(_h[0], 2) * _y[2] - _h[1] * _y[0] * (2 * _h[0] + _h[1]) +
            _y[1] * std::pow(_h[0] + _h[1], 2)) /
           (_h[0] * _h[1] * (_h[0] + _h[1]));
 
  Real hlast = _h[_n_intervals - 1];
  Real hsecond = _h[_n_intervals - 2];
  Real ylast = _y[_n_knots - 1];
  Real ysecond = _y[_n_knots - 2];
  Real ythird = _y[_n_knots - 3];
  _yp[_n_knots - 1] = (hsecond * ylast * (hsecond + 2 * hlast) + std::pow(hlast, 2) * ythird -
                       ysecond * std::pow(hsecond + hlast, 2)) /
                      (hsecond * hlast * (hsecond + hlast));
}

Referenced by solve().

modify_derivs()

void MonotoneCubicInterpolation::modify_derivs ( const Real &  alpha, const Real &  beta, const Real &  delta, Real &  yp_lo, Real &  yp_hi  )

protectedvirtual

Definition at line 271 of file MonotoneCubicInterpolation.C.

{
  Real tau = 3. / std::sqrt(std::pow(alpha, 2) + std::pow(beta, 2));
  Real alpha_star = alpha * tau;
  Real beta_star = beta * tau;
  yp_lo = alpha_star * delta;
  yp_hi = beta_star * delta;
}

Referenced by solve().

p()

Real MonotoneCubicInterpolation::p ( const Real &  xhi, const Real &  xlo, const Real &  fhi, const Real &  flo, const Real &  dhi, const Real &  dlo, const Real &  x  ) const

protectedvirtual

Definition at line 210 of file MonotoneCubicInterpolation.C.

{
  return flo * h1(xhi, xlo, x) + fhi * h2(xhi, xlo, x) + dlo * h3(xhi, xlo, x) +
         dhi * h4(xhi, xlo, x);
}

Referenced by sample().

pDoublePrime()

Real MonotoneCubicInterpolation::pDoublePrime ( const Real &  xhi, const Real &  xlo, const Real &  fhi, const Real &  flo, const Real &  dhi, const Real &  dlo, const Real &  x  ) const

protectedvirtual

Definition at line 236 of file MonotoneCubicInterpolation.C.

{
  return flo * h1DoublePrime(xhi, xlo, x) + fhi * h2DoublePrime(xhi, xlo, x) +
         dlo * h3DoublePrime(xhi, xlo, x) + dhi * h4DoublePrime(xhi, xlo, x);
}

Referenced by sample2ndDerivative().

phi()

Real MonotoneCubicInterpolation::phi ( const Real &  t ) const

protected

Definition at line 70 of file MonotoneCubicInterpolation.C.

{
  return 3. * t * t - 2. * t * t * t;
}

Referenced by h1(), and h2().

phiDoublePrime()

Real MonotoneCubicInterpolation::phiDoublePrime ( const Real &  t ) const

protected

Definition at line 82 of file MonotoneCubicInterpolation.C.

{
  return 6. - 12. * t;
}

Referenced by h1DoublePrime(), and h2DoublePrime().

phiPrime()

Real MonotoneCubicInterpolation::phiPrime ( const Real &  t ) const

protected

Definition at line 76 of file MonotoneCubicInterpolation.C.

{
  return 6. * t - 6. * t * t;
}

Referenced by h1Prime(), and h2Prime().

pPrime()

Real MonotoneCubicInterpolation::pPrime ( const Real &  xhi, const Real &  xlo, const Real &  fhi, const Real &  flo, const Real &  dhi, const Real &  dlo, const Real &  x  ) const

protectedvirtual

Definition at line 223 of file MonotoneCubicInterpolation.C.

{
  return flo * h1Prime(xhi, xlo, x) + fhi * h2Prime(xhi, xlo, x) + dlo * h3Prime(xhi, xlo, x) +
         dhi * h4Prime(xhi, xlo, x);
}

Referenced by sampleDerivative().

psi()

Real MonotoneCubicInterpolation::psi ( const Real &  t ) const

protected

Definition at line 88 of file MonotoneCubicInterpolation.C.

{
  return t * t * t - t * t;
}

Referenced by h3(), and h4().

psiDoublePrime()

Real MonotoneCubicInterpolation::psiDoublePrime ( const Real &  t ) const

protected

Definition at line 100 of file MonotoneCubicInterpolation.C.

{
  return 6. * t - 2.;
}

Referenced by h3DoublePrime(), and h4DoublePrime().

psiPrime()

Real MonotoneCubicInterpolation::psiPrime ( const Real &  t ) const

protected

Definition at line 94 of file MonotoneCubicInterpolation.C.

{
  return 3. * t * t - 2. * t;
}

Referenced by h3Prime(), and h4Prime().

sample()

Real MonotoneCubicInterpolation::sample ( const Real &  x ) const

virtual

This function will take an independent variable input and will return the dependent variable based on the generated fit.

Definition at line 347 of file MonotoneCubicInterpolation.C.

{
  // sanity check (empty MontoneCubicInterpolations are constructable
  // so we cannot put this into the errorCheck)
  assert(_x.size() > 0);
 
  unsigned int klo, khi;
  findInterval(x, klo, khi);
  return p(_x[khi], _x[klo], _y[khi], _y[klo], _yp[khi], _yp[klo], x);
}

Referenced by dumpCSV().

sample2ndDerivative()

Real MonotoneCubicInterpolation::sample2ndDerivative ( const Real &  x ) const

virtual

This function will take an independent variable input and will return the second derivative of the dependent variable with respect to the independent variable based on the generated fit.

Note that this can be discontinous at the knots.

Definition at line 367 of file MonotoneCubicInterpolation.C.

{
  unsigned int klo, khi;
  findInterval(x, klo, khi);
  return pDoublePrime(_x[khi], _x[klo], _y[khi], _y[klo], _yp[khi], _yp[klo], x);
}

Referenced by dumpCSV().

sampleDerivative()

Real MonotoneCubicInterpolation::sampleDerivative ( const Real &  x ) const

virtual

This function will take an independent variable input and will return the derivative of the dependent variable with respect to the independent variable based on the generated fit.

Definition at line 359 of file MonotoneCubicInterpolation.C.

{
  unsigned int klo, khi;
  findInterval(x, klo, khi);
  return pPrime(_x[khi], _x[klo], _y[khi], _y[klo], _yp[khi], _yp[klo], x);
}

Referenced by dumpCSV().

setData()

void MonotoneCubicInterpolation::setData ( const std::vector< Real > &  x, const std::vector< Real > &  y  )

virtual

Method generally used when MonotoneCubicInterpolation object was created using the empty constructor.

Takes two vectors of points for which to apply the fit. One should be of the independent variable while the other should be of the dependent variable. These values should have a one-to-one correspondence, e.g. the vectors must be of the same size.

Definition at line 31 of file MonotoneCubicInterpolation.C.

{
  _x = x;
  _y = y;
  errorCheck();
  solve();
}

sign()

Real MonotoneCubicInterpolation::sign ( const Real &  x ) const

protected

Definition at line 59 of file MonotoneCubicInterpolation.C.

{
  if (x < 0)
    return -1;
  else if (x > 0)
    return 1;
  else
    return 0;
}

Referenced by solve().

solve()

void MonotoneCubicInterpolation::solve ( )

protectedvirtual

Definition at line 282 of file MonotoneCubicInterpolation.C.

{
  _n_knots = _x.size(), _n_intervals = _x.size() - 1, _internal_knots = _x.size() - 2;
  _h.resize(_n_intervals);
  _yp.resize(_n_knots);
  _delta.resize(_n_intervals);
  _alpha.resize(_n_intervals);
  _beta.resize(_n_intervals);
 
  for (unsigned int i = 0; i < _n_intervals; ++i)
    _h[i] = _x[i + 1] - _x[i];
 
  initialize_derivs();
  for (unsigned int i = 0; i < _n_intervals; ++i)
    _delta[i] = (_y[i + 1] - _y[i]) / _h[i];
  if (sign(_delta[0]) != sign(_yp[0]))
    _yp[0] = 0;
  if (sign(_delta[_n_intervals - 1]) != sign(_yp[_n_knots - 1]))
    _yp[_n_knots - 1] = 0;
  for (unsigned int i = 0; i < _internal_knots; ++i)
    if (sign(_delta[i + 1]) == 0 || sign(_delta[i]) == 0 || sign(_delta[i + 1]) != sign(_delta[i]))
      _yp[1 + i] = 0;
 
  for (unsigned int i = 0; i < _n_intervals; ++i)
  {
    // Test for zero slope
    if (_yp[i] == 0)
      _alpha[i] = 0;
    else if (_delta[i] == 0)
      _alpha[i] = 4;
    else
      _alpha[i] = _yp[i] / _delta[i];
 
    // Test for zero slope
    if (_yp[i + 1] == 0)
      _beta[i] = 0;
    else if (_delta[i] == 0)
      _beta[i] = 4;
    else
      _beta[i] = _yp[i + 1] / _delta[i];
 
    // check if outside region of monotonicity
    if (std::pow(_alpha[i], 2) + std::pow(_beta[i], 2) > 9.)
      modify_derivs(_alpha[i], _beta[i], _delta[i], _yp[i], _yp[i + 1]);
  }
}

Referenced by MonotoneCubicInterpolation(), and setData().

Member Data Documentation

_alpha

std::vector<Real> MonotoneCubicInterpolation::_alpha

protected

Definition at line 146 of file MonotoneCubicInterpolation.h.

Referenced by solve().

_beta

std::vector<Real> MonotoneCubicInterpolation::_beta

protected

Definition at line 147 of file MonotoneCubicInterpolation.h.

Referenced by solve().

_delta

std::vector<Real> MonotoneCubicInterpolation::_delta

protected

Definition at line 145 of file MonotoneCubicInterpolation.h.

Referenced by solve().

_h

std::vector<Real> MonotoneCubicInterpolation::_h

protected

Definition at line 143 of file MonotoneCubicInterpolation.h.

Referenced by initialize_derivs(), and solve().

_internal_knots

unsigned int MonotoneCubicInterpolation::_internal_knots

protected

Definition at line 151 of file MonotoneCubicInterpolation.h.

Referenced by solve().

_n_intervals

unsigned int MonotoneCubicInterpolation::_n_intervals

protected

Definition at line 150 of file MonotoneCubicInterpolation.h.

Referenced by initialize_derivs(), and solve().

_n_knots

unsigned int MonotoneCubicInterpolation::_n_knots

protected

Definition at line 149 of file MonotoneCubicInterpolation.h.

Referenced by findInterval(), initialize_derivs(), and solve().

_x

std::vector<Real> MonotoneCubicInterpolation::_x

protected

Definition at line 141 of file MonotoneCubicInterpolation.h.

Referenced by errorCheck(), findInterval(), getSampleSize(), sample(), sample2ndDerivative(), sampleDerivative(), setData(), and solve().

_y

std::vector<Real> MonotoneCubicInterpolation::_y

protected

Definition at line 142 of file MonotoneCubicInterpolation.h.

Referenced by errorCheck(), initialize_derivs(), sample(), sample2ndDerivative(), sampleDerivative(), setData(), and solve().

_yp

std::vector<Real> MonotoneCubicInterpolation::_yp

protected

Definition at line 144 of file MonotoneCubicInterpolation.h.

Referenced by initialize_derivs(), sample(), sample2ndDerivative(), sampleDerivative(), and solve().

---

The documentation for this class was generated from the following files:

• include/utils/MonotoneCubicInterpolation.h
• src/utils/MonotoneCubicInterpolation.C