From: kanschat Date: Wed, 16 Sep 2009 18:55:58 +0000 (+0000) Subject: Lagrange polynomials of degree larger than 10 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=dc62e93f262870fcc3b10f37752bcea0c2cfb2e1;p=dealii-svn.git Lagrange polynomials of degree larger than 10 git-svn-id: https://svn.dealii.org/trunk@19466 0785d39b-7218-0410-832d-ea1e28bc413d --- diff --git a/deal.II/base/include/base/polynomial.h b/deal.II/base/include/base/polynomial.h index 330a6c411f..8eda982788 100644 --- a/deal.II/base/include/base/polynomial.h +++ b/deal.II/base/include/base/polynomial.h @@ -352,9 +352,10 @@ namespace Polynomials * constructor. */ static - std::vector + void compute_coefficients (const unsigned int n, - const unsigned int support_point); + const unsigned int support_point, + std::vector& a); }; /** diff --git a/deal.II/base/source/polynomial.cc b/deal.II/base/source/polynomial.cc index df38a18697..e7611a19b2 100644 --- a/deal.II/base/source/polynomial.cc +++ b/deal.II/base/source/polynomial.cc @@ -366,17 +366,47 @@ namespace Polynomials LagrangeEquidistant::LagrangeEquidistant (const unsigned int n, const unsigned int support_point) - : - Polynomial(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 two (2, 1.); + + for (unsigned int k=0;k<=n;++k) + { + if (k != support_point) + { + two[0] = -1.*k/n; + Polynomial factor(two); + factor.scale(1.*n/(support_point - k)); + (*this) *= factor; + } + } + } + } + - std::vector + void LagrangeEquidistant::compute_coefficients (const unsigned int n, - const unsigned int support_point) + const unsigned int support_point, + std::vector& a) { - std::vector a (n+1); Assert(support_point