]> https://gitweb.dealii.org/ - dealii-svn.git/commitdiff
Cleanup of computations of derivatives in Lagrange polynomial.
authorkronbichler <kronbichler@0785d39b-7218-0410-832d-ea1e28bc413d>
Wed, 21 Sep 2011 12:09:28 +0000 (12:09 +0000)
committerkronbichler <kronbichler@0785d39b-7218-0410-832d-ea1e28bc413d>
Wed, 21 Sep 2011 12:09:28 +0000 (12:09 +0000)
git-svn-id: https://svn.dealii.org/trunk@24350 0785d39b-7218-0410-832d-ea1e28bc413d

deal.II/source/base/polynomial.cc

index 83f964bccc0d35db096da1f229b5d6acb19da33f..b15144ff82ed079a01ddbd5d9c6ba262bea2f4f0 100644 (file)
@@ -114,13 +114,13 @@ namespace Polynomials
                                // (x-x_1)*(x-x_2)*...*(x-x_n), expand the
                                // derivatives like automatic differentiation
                                // does.
-       values[0] = 1.;
-       for (unsigned int d=1; d<values_size; ++d)
-         values[d] = 0.;
        const unsigned int n_supp = lagrange_support_points.size();
        switch (values_size)
          {
          default:
+           values[0] = 1;
+           for (unsigned int d=1; d<values_size; ++d)
+             values[d] = 0;
            for (unsigned int i=0; i<n_supp; ++i)
              {
                const number v = x-lagrange_support_points[i];
@@ -131,13 +131,33 @@ namespace Polynomials
                                // product rule for the old value and the new
                                // variable 'v', i.e., expand value v and
                                // derivative one). since we reuse a value
-                               // from the next lower derivative, need to
-                               // start from the highest derivative
-               for (unsigned int d=values_size-1; d>0; --d)
-                 values[d] = (values[d] * v +
-                              static_cast<number>(d) * values[d-1]);
+                               // from the next lower derivative from the
+                               // steps before, need to start from the
+                               // highest derivative
+               for (unsigned int k=values_size-1; k>0; --k)
+                 values[k] = (values[k] * v + values[k-1]);
                values[0] *= v;
              }
+                               // finally, multiply by the weight in the
+                               // Lagrange denominator. Could be done instead
+                               // of setting values[0] = 1 above, but that
+                               // gives different accumulation of round-off
+                               // errors (multiplication is not associative)
+                               // compared to when we computed the weight,
+                               // and hence a basis function might not be
+                               // exactly one at the center point, which is
+                               // nice to have. We also multiply derivatives
+                               // by k! to transform the product p_n =
+                               // p^(n)(x)/k! into the actual form of the
+                               // derivative
+           {
+             number k_faculty = 1;
+             for (unsigned int k=0; k<values_size; ++k)
+               {
+                 values[k] *= k_faculty * lagrange_weight;
+                 k_faculty *= static_cast<number>(k+1);
+               }
+           }
            break;
 
                                // manually implement size 1 (values only),
@@ -146,43 +166,44 @@ namespace Polynomials
                                // might be called often. then, we can unroll
                                // the loop.
          case 1:
+           values[0] = 1;
            for (unsigned int i=0; i<n_supp; ++i)
              {
                const number v = x-lagrange_support_points[i];
                values[0] *= v;
              }
+           values[0] *= lagrange_weight;
            break;
+
          case 2:
+           values[0] = 1;
+           values[1] = 0;
            for (unsigned int i=0; i<n_supp; ++i)
              {
                const number v = x-lagrange_support_points[i];
                values[1] = values[1] * v + values[0];
                values[0] *= v;
              }
+           values[0] *= lagrange_weight;
+           values[1] *= lagrange_weight;
            break;
+
          case 3:
+           values[0] = 1;
+           values[1] = 0;
+           values[2] = 0;
            for (unsigned int i=0; i<n_supp; ++i)
              {
                const number v = x-lagrange_support_points[i];
-               values[2] = values[2] * v + static_cast<number>(2) * values[1];
+               values[2] = values[2] * v + values[1];
                values[1] = values[1] * v + values[0];
                values[0] *= v;
              }
+           values[0] *= lagrange_weight;
+           values[1] *= lagrange_weight;
+           values[2] *= static_cast<number>(2) * lagrange_weight;
            break;
          }
-
-                               // finally, multiply by the weight in the
-                               // Lagrange denominator. Could be done instead
-                               // of setting values[0] = 1 above, but that
-                               // gives different accumulation of round-off
-                               // errors (multiplication is not associative)
-                               // compared to when we computed the weight,
-                               // and hence a basis function might not be
-                               // exactly one at the center point, which is
-                               // nice to have
-       for (unsigned int d=0; d<values_size; ++d)
-         values[d] *= lagrange_weight;
-
        return;
       }
 

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.