]> https://gitweb.dealii.org/ - dealii-svn.git/commitdiff
polynomial multiplication and Lagrange class started
authorguido <guido@0785d39b-7218-0410-832d-ea1e28bc413d>
Mon, 30 May 2005 13:17:12 +0000 (13:17 +0000)
committerguido <guido@0785d39b-7218-0410-832d-ea1e28bc413d>
Mon, 30 May 2005 13:17:12 +0000 (13:17 +0000)
git-svn-id: https://svn.dealii.org/trunk@10776 0785d39b-7218-0410-832d-ea1e28bc413d

deal.II/base/include/base/polynomial.h
deal.II/base/source/polynomial.cc
tests/base/polynomial1d.cc
tests/results/i686-pc-linux-gnu+gcc3.2/base/polynomial1d.output

index 6c1365a1a2cf94b8b41c2c0bd69b2ff95d6aaa3b..9096b55076721fc7e104a4ff52a92cfc3c87f057 100644 (file)
@@ -21,6 +21,7 @@
 
 #include <vector>
 
+template <int dim> class Point;
 
 /**
  * A namespace in which classes relating to the description of
@@ -163,6 +164,11 @@ namespace Polynomials
                                        */
       Polynomial<number>& operator *= (const double s);
 
+                                      /**
+                                       * Multiply with another polynomial.
+                                       */
+      Polynomial<number>& operator *= (const Polynomial<number>& p);
+
                                       /**
                                        * Add a second polynomial.
                                        */
@@ -333,7 +339,31 @@ namespace Polynomials
                             const unsigned int support_point);
   };
 
-
+/**
+ * Lagrange polynomials for an arbistrary set of interpolation points.
+ *
+ * @author Guido Kanschat, 2005
+ */
+  class Lagrange
+  {
+    public:
+                                      /**
+                                       * Given a set of points, this
+                                       * function returns all
+                                       * Lagrange polynomials for
+                                       * interpolation of these
+                                       * points. The number of
+                                       * polynomials is equal to the
+                                       * number of points and the
+                                       * maximum degree is one less.
+                                       */
+      static
+      std::vector<Polynomial<double> >
+      generate_complete_basis (const std::vector<Point<1> >& points);
+  };
+  
+  
+  
 /**
  * Legendre polynomials of arbitrary degree on <tt>[0,1]</tt>.
  *
index a9d3f3faba1db7ee37093dae5576d07835a62047..1e3e6589bded07df0afb14781066ed3e4c7612e2 100644 (file)
@@ -13,6 +13,7 @@
 
 
 #include <base/polynomial.h>
+#include <base/point.h>
 #include <base/exceptions.h>
 #include <base/thread_management.h>
 
@@ -167,6 +168,24 @@ namespace Polynomials
   }
 
   
+  template <typename number>
+  Polynomial<number>&
+  Polynomial<number>::operator *= (const Polynomial<number>& p)
+  {
+                                    // Degree of the product
+    unsigned int new_degree = this->degree() + p.degree();
+
+    std::vector<number> new_coefficients(new_degree+1, 0.);
+    
+    for (unsigned int i=0; i<p.coefficients.size(); ++i)
+      for (unsigned int j=0; j<this->coefficients.size(); ++j)
+      new_coefficients[i+j] += this->coefficients[j]*p.coefficients[i];
+    this->coefficients = new_coefficients;
+    
+    return *this;
+  }
+
+  
   template <typename number>
   Polynomial<number>&
   Polynomial<number>::operator += (const Polynomial<number>& p)
@@ -589,6 +608,38 @@ namespace Polynomials
   }
 
 
+//----------------------------------------------------------------------//
+
+  
+  std::vector<Polynomial<double> >
+  Lagrange::generate_complete_basis (const std::vector<Point<1> >& points)
+  {
+    std::vector<Polynomial<double> > p(points.size());
+                                    // polynomials are built as
+                                    // products of linear
+                                    // factors. The coefficient in
+                                    // front of the linear term is
+                                    // always 1.
+    std::vector<double> linear(2, 1.);
+                                    // We start with a constant polynomial
+    std::vector<double> one(1, 1.);
+    
+    for (unsigned int i=0;i<p.size();++i)
+      {
+                                        // Construct interpolation formula
+       p[i] = Polynomial<double>(one);
+       for (unsigned int k=0;k<p.size();++k)
+         if (k != i)
+           {
+             linear[0] = -points[k](0);
+             Polynomial<double> factor(linear);
+             factor *= 1./(points[i](0)-points[k](0));
+             p[i] *= factor;
+           }
+      }
+    return p;
+  }
+  
 
 // ------------------ class Legendre --------------- //
 
index dcd7647e8f3faf944106ba9738b3b4ad7734f297..e2a5f782508c551065ad3e46f7404928088c1d4e 100644 (file)
@@ -77,6 +77,9 @@ void polynomial_arithmetic ()
   p1 *= 2.;
   std::cerr << "*2" << std::endl;
   p1.print(std::cerr);
+  std::cerr << "*P2" << std::endl;
+  p2 *= p1;
+  p2.print(std::cerr);
   
   for (unsigned int i=0;i<7;++i)
     {
index d73e32dc27e77429a4027f0808d6aaf5b287d12f..fb8dcb3a55afa96dbf6e1ec75b49f34216c7e4d3 100644 (file)
@@ -32,6 +32,16 @@ P1+P2+x^5
 7.40000 x^2
 7.40000 x^1
 4.80000 x^0
+*P2
+8.00000 x^8
+4.80000 x^7
+40.8000 x^6
+54.4000 x^5
+82.5600 x^4
+91.9200 x^3
+64.8000 x^2
+41.8400 x^1
+13.4400 x^0
 derive
 10.0000 x^4
 0.00000 x^3

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.