]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Initial implementation of new classes Polynomial (includes Horner scheme) and Lagrang...
authorRalf Hartmann <Ralf.Hartmann@dlr.de>
Fri, 15 Dec 2000 06:59:08 +0000 (06:59 +0000)
committerRalf Hartmann <Ralf.Hartmann@dlr.de>
Fri, 15 Dec 2000 06:59:08 +0000 (06:59 +0000)
git-svn-id: https://svn.dealii.org/trunk@3534 0785d39b-7218-0410-832d-ea1e28bc413d

deal.II/base/include/base/polynomial.h [new file with mode: 0644]
deal.II/base/source/polynomial.cc [new file with mode: 0644]

diff --git a/deal.II/base/include/base/polynomial.h b/deal.II/base/include/base/polynomial.h
new file mode 100644 (file)
index 0000000..2bc65b7
--- /dev/null
@@ -0,0 +1,95 @@
+/*----------------------------   polynomial.h     ---------------------------*/
+/*      $Id$                 */
+/*                Ralf Hartmann, University of Heidelberg                */
+#ifndef __polynomial_H
+#define __polynomial_H
+/*----------------------------   polynomial.h     ---------------------------*/
+
+
+#include <base/exceptions.h>
+
+#include <vector.h>
+
+/**
+ * Base class for all 1D polynomials.
+ *
+ * @author Ralf Hartmann, 2000
+ */
+class Polynomial
+{
+  public:
+                                    /**
+                                     * Constructor.
+                                     */
+    Polynomial(const vector<double> &a);
+
+                                    /**
+                                     * Returns the values and the
+                                     * derivatives of the @p{Polynomial}
+                                     * at point @p{x}. @p{values[i],
+                                     * i=0,...,values.size()}
+                                     * includes the @p{i}th
+                                     * derivative.
+                                     *
+                                     * This function uses the Horner
+                                     * scheme.
+                                     */
+    void value(double x, vector<double> &values) const;
+
+  protected:
+
+                                    /**
+                                     * Coefficients of the polynomial
+                                     * $\sum_ia_ix^i$. This vector is
+                                     * filled by the constructor of
+                                     * derived classes.
+                                     */
+    const vector<double> coefficients;
+};
+
+
+
+/**
+ * Class of Lagrange polynomials with equidistant interpolation
+ * points. The polynomial of order @p{n} has got @p{n+1} interpolation
+ * points. The interpolation points are x=0, x=1 and x=intermediate
+ * points in ]0,1[ in ascending order. This order gives an index to
+ * each interpolation point.  A Lagrangian polynomial equals 1 at one
+ * interpolation point that is called `support point', and 0 at all other
+ * interpolation points.
+ *
+ * @author Ralf Hartmann, 2000
+ */
+class LagrangeEquidistant: public Polynomial
+{
+  public:
+                                    /**
+                                     * Constructor. Takes the order
+                                     * @p{n} of the Lagrangian
+                                     * polynom and the index
+                                     * @p{support_point} of the
+                                     * support point. Fills the
+                                     * @p{coefficients} of the base
+                                     * class @p{Polynomial}.
+                                     */
+    LagrangeEquidistant(unsigned int n, unsigned int support_point);
+
+  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.
+                                     */
+    static vector<double> compute_coefficients(unsigned int n, unsigned int support_point);
+};
+
+
+
+/*----------------------------   polynomial.h     ---------------------------*/
+/* end of #ifndef __polynomial_H */
+#endif
+/*----------------------------   polynomial.h     ---------------------------*/
diff --git a/deal.II/base/source/polynomial.cc b/deal.II/base/source/polynomial.cc
new file mode 100644 (file)
index 0000000..62401eb
--- /dev/null
@@ -0,0 +1,180 @@
+//----------------------------  polynomial.cc  ---------------------------
+//      $Id$                 */
+//    Version: $Name$
+//
+//    Copyright (C) 2000 by the deal.II authors
+//
+//    This file is subject to QPL and may not be  distributed
+//    without copyright and license information. Please refer
+//    to the file deal.II/doc/license.html for the  text  and
+//    further information on this license.
+//
+//----------------------------   polynomial.cc     -----------------------
+
+
+#include<base/polynomial.h>
+
+
+Polynomial::Polynomial(const vector<double> &a):
+               coefficients(a)
+{}
+
+
+void Polynomial::value(double x, vector<double> &values) const
+{
+  const unsigned int m=coefficients.size();
+  vector<double> a(coefficients);
+                                  // Horner scheme
+  unsigned int j_faculty=1;
+  for (unsigned int j=0; j<values.size(); ++j)
+    {      
+      for (int k=m-1; k>=static_cast<int>(j); --k)
+       a[k]+=x*a[k+1];
+      values[j]=j_faculty*a[j];
+
+      j_faculty*=j+1;
+    }
+}
+
+
+// ------------------------------------------------------------ //
+
+
+
+LagrangeEquidistant::LagrangeEquidistant(unsigned int n, unsigned int support_point):
+               Polynomial(compute_coefficients(n,support_point))
+{}
+
+
+vector<double> LagrangeEquidistant::compute_coefficients(unsigned int n, unsigned int support_point)
+{
+  vector<double> a;
+  a.resize(n+1);
+  Assert(support_point<n+1, ExcIndexRange(support_point, 0, n+1));
+  switch (n)
+    {
+      case 0:
+           switch (support_point)
+             {
+               case 0:
+                     a[0]=1.;
+                     break;
+               default:
+                     Assert(false, ExcInternalError());
+             }
+           break;
+      case 1:
+           switch (support_point)
+             {
+               case 0:
+                     a[0]=1.;
+                     a[1]=-1.;
+                     break;
+               case 1:
+                     a[0]=0.;
+                     a[1]=1.;
+                     break;
+               default:
+                     Assert(false, ExcInternalError());
+             }
+           break;
+      case 2:
+           switch (support_point)
+             {
+               case 0:
+                     a[0]=1.;
+                     a[1]=-3.;
+                     a[2]=2.;
+                     break;
+               case 1:
+                     a[0]=0.;
+                     a[1]=-1.;               
+                     a[2]=2.;
+                     break;
+               case 2:
+                     a[0]=0.;
+                     a[1]=4.;                
+                     a[2]=-4.;
+                     break;
+               default:
+                     Assert(false, ExcInternalError());
+             }
+           break;
+      case 3:
+           switch (support_point)
+             {
+               case 0:
+                     a[0]=1.0;
+                     a[1]=-11.0/2.0;
+                     a[2]=9.0;
+                     a[3]=-9.0/2.0;
+                     break;
+               case 1:
+                     a[0]=0.;
+                     a[1]=1.;                
+                     a[2]=-9.0/2.0;
+                     a[3]=9.0/2.0;
+                     break;
+               case 2:
+                     a[0]=0.;
+                     a[1]=9.0;               
+                     a[2]=-45.0/2.0;
+                     a[3]=27.0/2.0;
+                     break;
+               case 3:
+                     a[0]=0.;
+                     a[1]=-9.0/2.0;                  
+                     a[2]=18.0;
+                     a[3]=-27.0/2.0;
+                     break;
+               default:
+                     Assert(false, ExcInternalError());
+             }
+           break;
+      case 4:
+           switch (support_point)
+             {
+               case 0:
+                     a[0]=1.;
+                     a[1]=-25.0/3.0;
+                     a[2]=70.0/3.0;
+                     a[3]=-80.0/3.0;
+                     a[4]=32.0/3.0;
+                     break;
+               case 1:
+                     a[0]=0.;
+                     a[1]=-1.;
+                     a[2]=22.0/3.0;
+                     a[3]=-16.0;
+                     a[4]=32.0/3.0;
+                     break;
+               case 2:
+                     a[0]=0.;
+                     a[1]=16.0;
+                     a[2]=-208.0/3.0;
+                     a[3]=96.0;
+                     a[4]=-128.0/3.0;
+                     break;
+               case 3:
+                     a[0]=0.;
+                     a[1]=-12.0;
+                     a[2]=76.0;
+                     a[3]=-128.0;
+                     a[4]=64.0;
+                     break;
+               case 4:
+                     a[0]=0.;
+                     a[1]=16.0/3.0;
+                     a[2]=-112.0/3.0;
+                     a[3]=224.0/3.0;
+                     a[4]=-128.0/3.0;
+                     break;
+               default:
+                     Assert(false, ExcInternalError());
+             }
+           break;
+      default:
+           Assert(false, ExcNotImplemented());
+    }
+  return a;
+}

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.