#include <vector.h>
/**
- * Base class for all 1D polynomials.
+ * Base class for all 1D polynomials. A pollynomial is represented in
+ * this class by its coefficients, which are set through the
+ * constructor or by derived classes. Evaluation of a polynomial
+ * happens through the Horner scheme which provides both numerical
+ * stability and a minimal number of numerical operations.
*
* @author Ralf Hartmann, 2000
*/
{
public:
/**
- * Constructor.
+ * Constructor. The coefficients
+ * of the polynomial are passed
+ * as arguments, and denote the
+ * polynomial @p{\sum_i a[i]
+ * x^i}, i.e. the first element
+ * of the array denotes the
+ * constant term, the second the
+ * linear one, and so on. The
+ * order of the polynomial
+ * represented by this object is
+ * thus the number of elements in
+ * the @p{coefficient} array
+ * minus one.
*/
- Polynomial(const vector<double> &a);
+ Polynomial (const vector<double> &coefficients);
/**
- * Returns the values and the
- * derivatives of the @p{Polynomial}
- * at point @p{x}. @p{values[i],
+ * Return the value of this
+ * polynomial at the given point.
+ *
+ * This function uses the Horner
+ * scheme for numerical stability
+ * of the evaluation.
+ */
+ double value (const double x) const;
+
+ /**
+ * Return the values and the
+ * derivatives of the
+ * @p{Polynomial} at point @p{x}.
+ * @p{values[i],
* i=0,...,values.size()-1}
* includes the @p{i}th
- * derivative.
+ * derivative. The number of
+ * derivatives to be computed is
+ * thus determined by the size of
+ * the array passed.
*
* This function uses the Horner
- * scheme.
+ * scheme for numerical stability
+ * of the evaluation.
*/
- void value(double x, vector<double> &values) const;
+ void value (const double x,
+ vector<double> &values) const;
+ /**
+ * Exception
+ */
+ DeclException0 (ExcEmptyArray);
+
protected:
/**
* Coefficients of the polynomial
- * $\sum_ia_ix^i$. This vector is
- * filled by the constructor of
- * derived classes.
+ * $\sum_i a_i x^i$. This vector
+ * is filled by the constructor
+ * of this class and may be
+ * passed down by derived
+ * classes.
*/
const vector<double> coefficients;
};
* order. This order gives an index to each interpolation point. A
* Lagrangian polynomial equals 1 at one interpolation point that is
* then called `support point', and 0 at all other interpolation
- * points.
+ * points. For example, if the order is 3, and the support point is 1,
+ * then the polynomial represented by this object is of cubic and its
+ * value is 1 at the point @p{x=1/3}, and zero at the point @p{x=0},
+ * @p{x=2/3}, and @p{x=1}.
*
* @author Ralf Hartmann, 2000
*/
* @p{coefficients} of the base
* class @p{Polynomial}.
*/
- LagrangeEquidistant(unsigned int n, unsigned int support_point);
+ LagrangeEquidistant (const unsigned int n,
+ const unsigned int support_point);
+
+ /**
+ * Exception
+ */
+ DeclException1 (ExcInvalidSupportPoint,
+ int,
+ << "The support point " << arg1 << " is invalid.");
private:
* Computes the @p{coefficients}
* of the base class
* @p{Polynomial}. This function
- * is static to allow the
- * @p{coefficients} to be a
- * @p{const} vector.
+ * is @p{static} to allow to be
+ * called in the
+ * constructor. This in turn
+ * enables us to have the
+ * @p{coefficients} of the base
+ * class to be a @p{const}
+ * vector.
*/
- static vector<double> compute_coefficients(unsigned int n, unsigned int support_point);
+ static
+ vector<double>
+ compute_coefficients (const unsigned int n,
+ const unsigned int support_point);
};
#include <base/polynomial.h>
-Polynomial::Polynomial(const vector<double> &a):
+Polynomial::Polynomial (const vector<double> &a):
coefficients(a)
{}
-void Polynomial::value(double x, vector<double> &values) const
+
+double Polynomial::value (const double x) const
{
const unsigned int m=coefficients.size();
- vector<double> a(coefficients);
+
// Horner scheme
+ double value = coefficients.back();
+ for (int k=m-2; k>=0; --k)
+ value = value*x + coefficients[k];
+
+ return value;
+}
+
+
+
+void Polynomial::value (const double x,
+ vector<double> &values) const
+{
+ Assert (values.size() > 0, ExcEmptyArray());
+ const unsigned int m=coefficients.size();
+
+ // if we only need the value, then
+ // call the other function since
+ // that is significantly faster
+ // (there is no need to allocate
+ // and free memory, which is really
+ // expensive compared to all the
+ // other operations!)
+ if (m == 1)
+ {
+ values[0] = value(x);
+ return;
+ };
+
+ // if there are derivatives needed,
+ // then do it properly by the
+ // full Horner scheme
+ vector<double> a(coefficients);
unsigned int j_faculty=1;
for (unsigned int j=0; j<values.size(); ++j)
{
-LagrangeEquidistant::LagrangeEquidistant(unsigned int n, unsigned int support_point):
+LagrangeEquidistant::LagrangeEquidistant (const unsigned int n,
+ const unsigned int support_point):
Polynomial(compute_coefficients(n,support_point))
{}
-vector<double> LagrangeEquidistant::compute_coefficients(unsigned int n, unsigned int support_point)
+
+vector<double>
+LagrangeEquidistant::compute_coefficients (const unsigned int n,
+ const unsigned int support_point)
{
- vector<double> a;
- a.resize(n+1);
+ vector<double> a (n+1);
Assert(support_point<n+1, ExcIndexRange(support_point, 0, n+1));
+
switch (n)
{
case 0:
a[0]=1.;
break;
default:
- Assert(false, ExcInternalError());
+ Assert(false, ExcInvalidSupportPoint(support_point));
}
break;
case 1:
a[1]=1.;
break;
default:
- Assert(false, ExcInternalError());
+ Assert(false, ExcInvalidSupportPoint(support_point));
}
break;
case 2:
a[2]=2.;
break;
default:
- Assert(false, ExcInternalError());
+ Assert(false, ExcInvalidSupportPoint(support_point));
}
break;
case 3:
a[3]=9.0/2.0;
break;
default:
- Assert(false, ExcInternalError());
+ Assert(false, ExcInvalidSupportPoint(support_point));
}
break;
case 4:
a[4]=32.0/3.0;
break;
default:
- Assert(false, ExcInternalError());
+ Assert(false, ExcInvalidSupportPoint(support_point));
}
break;
default: