From 9770792e426e76b01d70778a45d0c9c25369b59a Mon Sep 17 00:00:00 2001 From: Wolfgang Bangerth Date: Sun, 17 Dec 2000 14:16:59 +0000 Subject: [PATCH] Add another value function to the polynomial class and test it. git-svn-id: https://svn.dealii.org/trunk@3544 0785d39b-7218-0410-832d-ea1e28bc413d --- deal.II/base/include/base/polynomial.h | 93 +++++++++++++++++++++----- deal.II/base/source/polynomial.cc | 61 +++++++++++++---- tests/base/polynomial_test.cc | 12 ++++ 3 files changed, 136 insertions(+), 30 deletions(-) diff --git a/deal.II/base/include/base/polynomial.h b/deal.II/base/include/base/polynomial.h index ca599d3deb..e5acb1986f 100644 --- a/deal.II/base/include/base/polynomial.h +++ b/deal.II/base/include/base/polynomial.h @@ -20,7 +20,11 @@ #include /** - * 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 */ @@ -28,30 +32,65 @@ class Polynomial { 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 &a); + Polynomial (const vector &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 &values) const; + void value (const double x, + vector &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 coefficients; }; @@ -65,7 +104,10 @@ class Polynomial * 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 */ @@ -81,7 +123,15 @@ class LagrangeEquidistant: public Polynomial * @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: @@ -89,11 +139,18 @@ class LagrangeEquidistant: public Polynomial * 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 compute_coefficients(unsigned int n, unsigned int support_point); + static + vector + compute_coefficients (const unsigned int n, + const unsigned int support_point); }; diff --git a/deal.II/base/source/polynomial.cc b/deal.II/base/source/polynomial.cc index 80ad3e5763..d55a8a8629 100644 --- a/deal.II/base/source/polynomial.cc +++ b/deal.II/base/source/polynomial.cc @@ -15,16 +15,49 @@ #include -Polynomial::Polynomial(const vector &a): +Polynomial::Polynomial (const vector &a): coefficients(a) {} -void Polynomial::value(double x, vector &values) const + +double Polynomial::value (const double x) const { const unsigned int m=coefficients.size(); - vector 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 &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 a(coefficients); unsigned int j_faculty=1; for (unsigned int j=0; j &values) const -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 LagrangeEquidistant::compute_coefficients(unsigned int n, unsigned int support_point) + +vector +LagrangeEquidistant::compute_coefficients (const unsigned int n, + const unsigned int support_point) { - vector a; - a.resize(n+1); + vector a (n+1); Assert(support_point LagrangeEquidistant::compute_coefficients(unsigned int n, unsigne a[0]=1.; break; default: - Assert(false, ExcInternalError()); + Assert(false, ExcInvalidSupportPoint(support_point)); } break; case 1: @@ -75,7 +112,7 @@ vector LagrangeEquidistant::compute_coefficients(unsigned int n, unsigne a[1]=1.; break; default: - Assert(false, ExcInternalError()); + Assert(false, ExcInvalidSupportPoint(support_point)); } break; case 2: @@ -97,7 +134,7 @@ vector LagrangeEquidistant::compute_coefficients(unsigned int n, unsigne a[2]=2.; break; default: - Assert(false, ExcInternalError()); + Assert(false, ExcInvalidSupportPoint(support_point)); } break; case 3: @@ -128,7 +165,7 @@ vector LagrangeEquidistant::compute_coefficients(unsigned int n, unsigne a[3]=9.0/2.0; break; default: - Assert(false, ExcInternalError()); + Assert(false, ExcInvalidSupportPoint(support_point)); } break; case 4: @@ -170,7 +207,7 @@ vector LagrangeEquidistant::compute_coefficients(unsigned int n, unsigne a[4]=32.0/3.0; break; default: - Assert(false, ExcInternalError()); + Assert(false, ExcInvalidSupportPoint(support_point)); } break; default: diff --git a/tests/base/polynomial_test.cc b/tests/base/polynomial_test.cc index bdc4267e74..b9bbe4261a 100644 --- a/tests/base/polynomial_test.cc +++ b/tests/base/polynomial_test.cc @@ -58,6 +58,18 @@ int main(int, char) else deallog << " false"; deallog << endl; + + // now also check + // whether the other + // @p{value} function + // returns the same + // result + if (polynom.value(x) != values[0]) + { + deallog << "The two `value' functions return different results!" + << endl; + abort (); + }; } } } -- 2.39.5