From 2e1143b593979a2b45594c2fd0eb14c22a4b2356 Mon Sep 17 00:00:00 2001 From: kronbichler Date: Wed, 21 Sep 2011 12:09:28 +0000 Subject: [PATCH] Cleanup of computations of derivatives in Lagrange polynomial. git-svn-id: https://svn.dealii.org/trunk@24350 0785d39b-7218-0410-832d-ea1e28bc413d --- deal.II/source/base/polynomial.cc | 65 ++++++++++++++++++++----------- 1 file changed, 43 insertions(+), 22 deletions(-) diff --git a/deal.II/source/base/polynomial.cc b/deal.II/source/base/polynomial.cc index 83f964bccc..b15144ff82 100644 --- a/deal.II/source/base/polynomial.cc +++ b/deal.II/source/base/polynomial.cc @@ -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; d0; --d) - values[d] = (values[d] * v + - static_cast(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(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(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(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