From b42f014d067d09d356e898fcfe2d6bd7aeec25d3 Mon Sep 17 00:00:00 2001
From: Andreas Lauser <and@poware.org>
Date: Tue, 21 Jun 2011 17:14:29 +0000
Subject: [PATCH] spline: add (optional) extrapolation in the eval and
eval-derivative methods
splines are extrapolated attaching straight lines at the spline's
endpoints which exhibits the same value and derivative as the spline
at the end points. this makes the eval() and evalDerivative() methods
more covenient to use.
git-svn-id: svn://svn.iws.uni-stuttgart.de/DUMUX/dumux/trunk@6053 2fb0f335-1f38-0410-981e-8018bf24f1b0
---
dumux/common/splinecommon_.hh | 144 ++++++++++++++++++++++------------
1 file changed, 95 insertions(+), 49 deletions(-)
diff --git a/dumux/common/splinecommon_.hh b/dumux/common/splinecommon_.hh
index 180cd5cfcc..608e90c3ee 100644
--- a/dumux/common/splinecommon_.hh
+++ b/dumux/common/splinecommon_.hh
@@ -128,65 +128,62 @@ public:
/*!
* \brief Evaluate the spline at a given position.
+ *
+ * \param x The value on the abscissa where the spline ought to be
+ * evaluated
+ * \param extrapolate If this parameter is set to true, the spline
+ * will be extended beyond its range by
+ * straight lines, if false calling extrapolate
+ * for \$f x \not [x_{min}, x_{max}]\f$ will
+ * cause a failed assertation.
*/
- Scalar eval(Scalar x) const
+ Scalar eval(Scalar x, bool extrapolate=false) const
{
- assert(applies(x));
-
- // See: J. Stoer: "Numerische Mathematik 1", 9th edition,
- // Springer, 2005, p. 109
-
- int i = segmentIdx_(x);
-
- Scalar h_i1 = h_(i + 1);
- Scalar x_i = x - x_(i);
- Scalar x_i1 = x_(i+1) - x;
-
- Scalar A_i =
- (y_(i+1) - y_(i))/h_i1
- -
- h_i1/6*(moment_(i+1) - moment_(i));
- Scalar B_i = y_(i) - moment_(i)* (h_i1*h_i1) / 6;
+ assert(extrapolate || applies(x));
+
+ // handle extrapolation
+ if (extrapolate) {
+ if (x < xMin()) {
+ Scalar m = evalDerivative(xMin(), /*segmentIdx=*/0);
+ Scalar y0 = y_(0);
+ return y0 + m*(x - xMin());
+ }
+ else if (x > xMax()) {
+ Scalar m = evalDerivative(xMax(), /*segmentIdx=*/numSamples_()-1);
+ Scalar y0 = y_(numSamples_()-1);
+ return y0 + m*(x - xMax());
+ }
+ }
- return
- moment_(i)* x_i1*x_i1*x_i1 / (6 * h_i1)
- +
- moment_(i + 1)* x_i*x_i*x_i / (6 * h_i1)
- +
- A_i*x_i
- +
- B_i;
+ return eval_(x, segmentIdx_(x));
}
/*!
* \brief Evaluate the spline's derivative at a given position.
+ *
+ * \param x The value on the abscissa where the spline's
+ * derivative ought to be evaluated
+ *
+ * \param extrapolate If this parameter is set to true, the spline
+ * will be extended beyond its range by
+ * straight lines, if false calling extrapolate
+ * for \$f x \not [x_{min}, x_{max}]\f$ will
+ * cause a failed assertation.
+
*/
- Scalar evalDerivative(Scalar x) const
+ Scalar evalDerivative(Scalar x, bool extrapolate=false) const
{
- assert(applies(x));
-
- // See: J. Stoer: "Numerische Mathematik 1", 9th edition,
- // Springer, 2005, p. 109
-
- int i = segmentIdx_(x);
-
- Scalar h_i1 = h_(i + 1);
- Scalar x_i = x - x_(i);
- Scalar x_i1 = x_(i+1) - x;
-
- Scalar A_i =
- (y_(i+1) - y_(i))/h_i1
- -
- h_i1/6*(moment_(i+1) - moment_(i));
-
- return
- -moment_(i) * x_i1*x_i1 / (2 * h_i1)
- +
- moment_(i + 1) * x_i*x_i / (2 * h_i1)
- +
- A_i;
+ assert(extrapolate || applies(x));
+ if (extrapolate) {
+ if (x < xMin())
+ evalDerivative_(xMin(), 0);
+ else if (x > xMax())
+ evalDerivative_(xMax(), numSamples_() - 1);
+ }
+
+ return evalDerivative_(x, segmentIdx_(x));
}
-
+
/*!
* \brief Find the intersections of the spline with a cubic
* polynomial in the whole intervall, throws
@@ -484,6 +481,55 @@ protected:
d[n] = 0.0;
}
+ // evaluate the spline at a given the position and given the
+ // segment index
+ Scalar eval_(Scalar x, int i) const
+ {
+ // See: J. Stoer: "Numerische Mathematik 1", 9th edition,
+ // Springer, 2005, p. 109
+ Scalar h_i1 = h_(i + 1);
+ Scalar x_i = x - x_(i);
+ Scalar x_i1 = x_(i+1) - x;
+
+ Scalar A_i =
+ (y_(i+1) - y_(i))/h_i1
+ -
+ h_i1/6*(moment_(i+1) - moment_(i));
+ Scalar B_i = y_(i) - moment_(i)* (h_i1*h_i1) / 6;
+
+ return
+ moment_(i)* x_i1*x_i1*x_i1 / (6 * h_i1)
+ +
+ moment_(i + 1)* x_i*x_i*x_i / (6 * h_i1)
+ +
+ A_i*x_i
+ +
+ B_i;
+ }
+
+ // evaluate the derivative of a spline given the actual position
+ // and the segment index
+ Scalar evalDerivative_(Scalar x, int i) const
+ {
+ // See: J. Stoer: "Numerische Mathematik 1", 9th edition,
+ // Springer, 2005, p. 109
+ Scalar h_i1 = h_(i + 1);
+ Scalar x_i = x - x_(i);
+ Scalar x_i1 = x_(i+1) - x;
+
+ Scalar A_i =
+ (y_(i+1) - y_(i))/h_i1
+ -
+ h_i1/6*(moment_(i+1) - moment_(i));
+
+ return
+ -moment_(i) * x_i1*x_i1 / (2 * h_i1)
+ +
+ moment_(i + 1) * x_i*x_i / (2 * h_i1)
+ +
+ A_i;
+ }
+
// returns the monotonicality of an interval of a spline segment
//
--
GitLab