]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Lagrange polynomials of degree larger than 10
authorGuido Kanschat <dr.guido.kanschat@gmail.com>
Wed, 16 Sep 2009 18:55:58 +0000 (18:55 +0000)
committerGuido Kanschat <dr.guido.kanschat@gmail.com>
Wed, 16 Sep 2009 18:55:58 +0000 (18:55 +0000)
git-svn-id: https://svn.dealii.org/trunk@19466 0785d39b-7218-0410-832d-ea1e28bc413d

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

index 330a6c411fcc4010a4f00ea002d9172392e43814..8eda9827884d14780e4a5546dda5d5733b6fb7fd 100644 (file)
@@ -352,9 +352,10 @@ namespace Polynomials
                                         * constructor.
                                         */
       static 
-      std::vector<double> 
+      void
       compute_coefficients (const unsigned int n,
-                            const unsigned int support_point);
+                            const unsigned int support_point,
+                           std::vector<double>& a);
   };
 
 /**
index df38a18697f194c651238dcf6e6db64d4b83f7c2..e7611a19b290c81a75cfc629d3b2c15be5d65eae 100644 (file)
@@ -366,17 +366,47 @@ namespace Polynomials
 
   LagrangeEquidistant::LagrangeEquidistant (const unsigned int n,
                                             const unsigned int support_point)
-                 :
-                  Polynomial<double>(compute_coefficients(n,support_point))
-  {}
-
-
+  {
+    if (n <= 10)
+      {
+       this->coefficients.resize(n+1);
+       compute_coefficients(n, support_point, this->coefficients);
+      }
+    else
+      {
+                                        // We have precomputed tables
+                                        // up to degree 10. For
+                                        // higher order, we have to
+                                        // compute by hand.
+       
+                                        // Start with the constant one
+       this->coefficients.resize(1);
+       this->coefficients[0] = 1.;
+
+                                        // Then compute the Lagrange
+                                        // polynomial as the product
+                                        // of linear factors
+       std::vector<double> two (2, 1.);
+       
+       for (unsigned int k=0;k<=n;++k)
+         {
+           if (k != support_point)
+             {
+               two[0] = -1.*k/n;
+               Polynomial<double> factor(two);
+               factor.scale(1.*n/(support_point - k));
+               (*this) *= factor;
+             }
+         }  
+      }
+  }
+  
 
-  std::vector<double> 
+  void
   LagrangeEquidistant::compute_coefficients (const unsigned int n,
-                                             const unsigned int support_point)
+                                             const unsigned int support_point,
+                                            std::vector<double>& a)
   {
-    std::vector<double> a (n+1);
     Assert(support_point<n+1, ExcIndexRange(support_point, 0, n+1));
 
     unsigned int n_functions=n+1;
@@ -583,14 +613,12 @@ namespace Polynomials
           break;
         }
         default:
-              Assert(false, ExcNotImplemented());
+             Assert(false, ExcInternalError())
       }
 
     Assert(x!=0, ExcInternalError());
     for (unsigned int i=0; i<n_functions; ++i)
       a[i]=x[support_point*n_functions+i];
-  
-    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.