]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Add another value function to the polynomial class and test it.
authorWolfgang Bangerth <bangerth@math.tamu.edu>
Sun, 17 Dec 2000 14:16:59 +0000 (14:16 +0000)
committerWolfgang Bangerth <bangerth@math.tamu.edu>
Sun, 17 Dec 2000 14:16:59 +0000 (14:16 +0000)
git-svn-id: https://svn.dealii.org/trunk@3544 0785d39b-7218-0410-832d-ea1e28bc413d

deal.II/base/include/base/polynomial.h
deal.II/base/source/polynomial.cc
tests/base/polynomial_test.cc

index ca599d3debd0b48950939603433434bf4956d2a7..e5acb1986f272c7e8b58558e13a71f94550c1788 100644 (file)
 #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
  */
@@ -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<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;
 };
@@ -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<double> compute_coefficients(unsigned int n, unsigned int support_point);
+    static 
+    vector<double> 
+    compute_coefficients (const unsigned int n,
+                         const unsigned int support_point);
 };
 
 
index 80ad3e57632e5e1d5c459ce4d7138d9574ea95f4..d55a8a862939d6947ebb513b351a040409669cad 100644 (file)
 #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)
     {      
@@ -41,16 +74,20 @@ void Polynomial::value(double x, vector<double> &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<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:
@@ -60,7 +97,7 @@ vector<double> 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<double> 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<double> 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<double> 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<double> 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:
index bdc4267e748f107f9f59d855d2e189e76cc2af06..b9bbe4261a3e1792486718ee84e5daccf21f844d 100644 (file)
@@ -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 ();
+               };
            }
        }
     }

In the beginning the Universe was created. This has made a lot of people very angry and has been widely regarded as a bad move.

Douglas Adams


Typeset in Trocchi and Trocchi Bold Sans Serif.