]> https://gitweb.dealii.org/ - dealii-svn.git/commitdiff
Implemented Legendre-Gauss-Lobatto quadrature.
authorprill <prill@0785d39b-7218-0410-832d-ea1e28bc413d>
Thu, 2 Nov 2006 14:40:04 +0000 (14:40 +0000)
committerprill <prill@0785d39b-7218-0410-832d-ea1e28bc413d>
Thu, 2 Nov 2006 14:40:04 +0000 (14:40 +0000)
git-svn-id: https://svn.dealii.org/trunk@14140 0785d39b-7218-0410-832d-ea1e28bc413d

deal.II/base/include/base/quadrature_lib.h
deal.II/base/source/quadrature_lib.cc
deal.II/doc/news/changes.html

index 6aa72a8e124230b151c44e1c10b0e7774b0c3b65..c34825041df5f9f31bf8abde2fb2e2f642ab0513 100644 (file)
@@ -49,19 +49,27 @@ class QGauss : public Quadrature<dim>
  * The Gauss-Lobatto quadrature rule.
  *
  * This modification of the Gauss quadrature uses the two interval end
- * points as well. Being exact for polynomials of degree <i>2n-2</i>,
- * this formula is suboptimal by one degree.
+ * points as well. Being exact for polynomials of degree <i>2n-3</i>,
+ * this formula is suboptimal by two degrees.
  *
- * The quadrature points are interval end points plus the zeroes of
+ * The quadrature points are interval end points plus the roots of
  * the derivative of the Legendre polynomial <i>P<sub>n-1</sub></i> of
  * degree <i>n-1</i>. The quadrature weights are
  * <i>2/(n(n-1)(P<sub>n-1</sub>(x<sub>i</sub>)<sup>2</sup>)</i>.
  *
- * @note The quadrature weights are not implemented yet.
+ * Note: This implementation has not yet been optimized concerning
+ *       numerical stability and efficiency. It can be easily adapted
+ *       to the general case of Gauss-Lobatto-Jacobi-Bouzitat quadrature
+ *       with arbitrary parameters <i>alpha</i>, <i>beta</i>, of which
+ *       the Gauss-Lobatto-Legendre quadrature (<i>alpha = beta = 0</i>)
+ *       is a special case.
  *
- * @sa http://en.wikipedia.org/wiki/Handbook_of_Mathematical_Functions
+ * @sa http://en.wikipedia.org/wiki/Handbook_of_Mathematical_Functions 
+ * @sa Karniadakis, G.E. and Sherwin, S.J.:
+ *     Spectral/hp element methods for computational fluid dynamics. 
+ *     Oxford: Oxford University Press, 2005 
  *
- * @author Guido Kanschat, 2005, 2006
+ * @author Guido Kanschat, 2005, 2006; F. Prill, 2006
  */
 template<int dim>
 class QGaussLobatto : public Quadrature<dim>
@@ -73,6 +81,47 @@ class QGaussLobatto : public Quadrature<dim>
                                      * (in each space direction).
                                      */
     QGaussLobatto(const unsigned int n);
+
+  protected:
+                                    /**
+                                     * Compute Legendre-Gauss-Lobatto quadrature
+                                     * points in the interval [-1, +1].
+                                     * @param q  number of points.
+                                     * @return vector containing nodes.
+                                     */
+    std::vector<double> compute_quadrature_points(const unsigned int q,
+                                                  const int alpha,
+                                                  const int beta) const;
+
+                                    /**
+                                     * Compute Legendre-Gauss-Lobatto quadrature
+                                     * weights.
+                                     * @param x  quadrature points.
+                                     * @return vector containing weights.
+                                     */
+    std::vector<double> compute_quadrature_weights(std::vector<double>& x,
+                                                   const int alpha,
+                                                   const int beta) const;
+    
+                                    /**
+                                     * Evaluate a Jacobi polynomial
+                                     * \f$ P^{\alpha, \beta}_n(x) \f$. 
+                                     * Note: The Jacobi polynomials are
+                                     * not orthonormal and defined on
+                                     * the interval [-1, +1].
+                                     * @param x  point of evaluation.
+                                     */
+    double JacobiP(const double x,
+                   const int alpha,
+                   const int beta,
+                   const unsigned int n) const;
+
+                                    /**
+                                     * Evaluate the Gamma function
+                                     * \f[ \Gamma(n) = (n-1)! \f]. 
+                                     * @param n  point of evaluation (integer).
+                                     */
+    unsigned int gamma(const unsigned int n) const;
 };
 
   
@@ -250,6 +299,18 @@ class QWeddle : public Quadrature<dim>
 
 template <> QGauss<1>::QGauss (const unsigned int n);
 template <> QGaussLobatto<1>::QGaussLobatto (const unsigned int n);
+template <>
+std::vector<double> QGaussLobatto<1>::
+compute_quadrature_points(const unsigned int, const int, const int) const;
+template <>
+std::vector<double> QGaussLobatto<1>::
+compute_quadrature_weights(std::vector<double>&, const int, const int) const;
+template <>
+double QGaussLobatto<1>::
+JacobiP(const double, const int, const int, const unsigned int) const;
+template <>
+unsigned int QGaussLobatto<1>::
+QGaussLobatto<1>::gamma(const unsigned int n) const;
 template <> QGauss2<1>::QGauss2 ();
 template <> QGauss3<1>::QGauss3 ();
 template <> QGauss4<1>::QGauss4 ();
index 223765e532a828ec6a71087365795742d656a7d0..49a41311dccffad18b0abdcb2f62b4031c89572c 100644 (file)
@@ -156,13 +156,139 @@ QGaussLobatto<1>::QGaussLobatto (unsigned int n)
                :
                Quadrature<1> (n)
 {
-  for (unsigned int k=0;k<n;++k)
+  std::vector<double> points  = compute_quadrature_points(n, 1, 1);
+  std::vector<double> w       = compute_quadrature_weights(points, 0, 0);
+
+                                  // scale points to the interval
+                                  // [0.0, 1.0]:
+  for (unsigned int i=0; i<points.size(); ++i)
+    {
+      this->quadrature_points[i] = Point<1>(0.5 + 0.5*points[i]);
+      this->weights[i]           = 0.5*w[i];
+    }
+}
+
+
+
+template <>
+std::vector<double> QGaussLobatto<1>::
+compute_quadrature_points(const unsigned int q,
+                         const int alpha,
+                         const int beta) const
+{
+  const unsigned int m = q-2; // no. of inner points
+  std::vector<double> x(m);
+  
+                                  // compute quadrature points with
+                                  // a Newton algorithm.
+  const double epsilon = 1.e-10; // tolerance
+  
+                                  // we take the zeros of the Chebyshev
+                                  // polynomial (alpha=beta=-0.5) as
+                                  // initial values:
+  for (unsigned int i=0; i<m; ++i)
+    x[i] = - std::cos( (double) (2*i+1)/(2*m) * M_PI );
+  
+  double r, s, J_x, f, delta;
+  for (unsigned int k=0; k<m; ++k)
+    {
+      r = x[k];
+      if (k>0)
+       r = (r + x[k-1])/2;
+
+      do 
+       {
+         s = 1.;
+         for (unsigned int i=0; i<k; ++i)
+           s /= r - x[i];
+
+         J_x   =  0.5*(alpha + beta + m + 1)*JacobiP(r, alpha+1, beta+1, m-1);
+         f     = JacobiP(r, alpha, beta, m);
+         delta = - f/(J_x - f*s);
+         r += delta;
+       }
+      while (std::fabs(delta) >= epsilon);
+      
+      x[k] = r;
+    } // for
+
+                                  // add boundary points:
+  x.insert(x.begin(), -1.);
+  x.push_back(+1.);
+
+  return x;
+}
+
+
+
+template <>
+std::vector<double> QGaussLobatto<1>::
+compute_quadrature_weights(std::vector<double>& x,
+                          const int alpha,
+                          const int beta) const
+{
+  const unsigned int q = x.size();
+  std::vector<double> w(q);
+  double s = 0;
+  
+  double factor = std::pow(2., alpha+beta+1)*gamma(alpha+q)*gamma(beta+q)/
+                 ((q-1)*gamma(q)*gamma(alpha+beta+q+1));
+  for (unsigned int i=0; i<q; ++i)
     {
-      this->quadrature_points[k] = Point<1>(.5-.5*std::cos(M_PI*k/(n-1)));
+      s = JacobiP(x[i], alpha, beta, q-1);
+      w[i] = factor/(s*s);
     }
+  w[0]   *= (beta + 1);
+  w[q-1] *= (alpha + 1);
+
+  return w;
+}
+
+
+
+template <>
+double QGaussLobatto<1>::JacobiP(const double x,
+                                const int alpha,
+                                const int beta,
+                                const unsigned int n) const
+{
+                                  // the Jacobi polynomial is evaluated
+                                  // using a recursion formula.
+  std::vector<double> p(n+1);
+  double v, a1, a2, a3, a4;
+  
+                                  // initial values P_0(x), P_1(x):
+  p[0] = 1.0;
+  if (n==0) return p[0];
+  p[1] = ((alpha+beta+2)*x + (alpha-beta))/2;
+  if (n==1) return p[1];
+  
+  for (unsigned int i=1; i<=(n-1); ++i)
+    {
+      v  = 2*i + alpha + beta;
+      a1 = 2*(i+1)*(i + alpha + beta + 1)*v;
+      a2 = (v + 1)*(alpha*alpha - beta*beta);
+      a3 = v*(v + 1)*(v + 2);
+      a4 = 2*(i+alpha)*(i+beta)*(v + 2);
+
+      p[i+1] =  1/a1*( (a2 + a3*x)*p[i] - a4*p[i-1]);
+    } // for
+  return p[n];
 }
 
 
+
+template <>
+unsigned int QGaussLobatto<1>::gamma(const unsigned int n) const
+{
+  unsigned int result = n - 1;
+  for (int i=n-2; i>1; --i)
+    result *= i;
+  return result;
+}
+
+
+
 template <>
 QGauss2<1>::QGauss2 ()
                 :
@@ -473,6 +599,13 @@ QGauss<dim>::QGauss (const unsigned int n)
 
 
 
+template <int dim>
+QGaussLobatto<dim>::QGaussLobatto (const unsigned int n)
+  :  Quadrature<dim> (QGaussLobatto<dim-1>(n), QGaussLobatto<1>(n))
+{}
+
+
+
 template <int dim>
 QGauss2<dim>::QGauss2 ()
                 :
index e041f93bc28e777151c08f26ac5fcc2fe53ed6fa..e44bf219deec7a51e74b235fc1a087fa7b873efb 100644 (file)
@@ -390,6 +390,11 @@ inconvenience this causes.
 <h3>base</h3>
 
 <ol>
+  <li> <p> Extended: The <code class="class">QGaussLobatto</code> quadrature rule computes
+           Legendre-Gauss-Lobatto nodes and quadrature weights.
+           <br>
+                (F. Prill 2006/11/02)
+       </p>
   <li> <p> Extended: The <code>contract</code> function family contracts two
            tensors. There is now a new version of this function, that contracts a
            tensor of rank three with a second one of rank two over given indices

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.