From 80f645959ebf02bad0cbc9075f023b9367a422c7 Mon Sep 17 00:00:00 2001 From: Ralf Hartmann Date: Fri, 15 Dec 2000 10:33:36 +0000 Subject: [PATCH] Lexicographical order of support points. git-svn-id: https://svn.dealii.org/trunk@3540 0785d39b-7218-0410-832d-ea1e28bc413d --- deal.II/base/include/base/polynomial.h | 12 +++---- deal.II/base/source/polynomial.cc | 40 +++++++++++------------ tests/base/polynomial_test.cc | 44 ++++++++------------------ tests/base/polynomial_test.checked | 38 +++++++--------------- 4 files changed, 52 insertions(+), 82 deletions(-) diff --git a/deal.II/base/include/base/polynomial.h b/deal.II/base/include/base/polynomial.h index 923f59db4c..ca599d3deb 100644 --- a/deal.II/base/include/base/polynomial.h +++ b/deal.II/base/include/base/polynomial.h @@ -36,7 +36,7 @@ class Polynomial * Returns the values and the * derivatives of the @p{Polynomial} * at point @p{x}. @p{values[i], - * i=0,...,values.size()} + * i=0,...,values.size()-1} * includes the @p{i}th * derivative. * @@ -61,11 +61,11 @@ class Polynomial /** * Class of Lagrange polynomials with equidistant interpolation * points. The polynomial of order @p{n} has got @p{n+1} interpolation - * points. The interpolation points are x=0, x=1 and x=intermediate - * points in ]0,1[ in ascending order. This order gives an index to - * each interpolation point. A Lagrangian polynomial equals 1 at one - * interpolation point that is called `support point', and 0 at all other - * interpolation points. + * points. The interpolation points are sorted in ascending + * order. This order gives an index to each interpolation point. A + * Lagrangian polynomial equals 1 at one interpolation point that is + * then called `support point', and 0 at all other interpolation + * points. * * @author Ralf Hartmann, 2000 */ diff --git a/deal.II/base/source/polynomial.cc b/deal.II/base/source/polynomial.cc index cacbecff75..80ad3e5763 100644 --- a/deal.II/base/source/polynomial.cc +++ b/deal.II/base/source/polynomial.cc @@ -88,13 +88,13 @@ vector LagrangeEquidistant::compute_coefficients(unsigned int n, unsigne break; case 1: a[0]=0.; - a[1]=-1.; - a[2]=2.; + a[1]=4.; + a[2]=-4.; break; case 2: a[0]=0.; - a[1]=4.; - a[2]=-4.; + a[1]=-1.; + a[2]=2.; break; default: Assert(false, ExcInternalError()); @@ -110,23 +110,23 @@ vector LagrangeEquidistant::compute_coefficients(unsigned int n, unsigne a[3]=-9.0/2.0; break; case 1: - a[0]=0.; - a[1]=1.; - a[2]=-9.0/2.0; - a[3]=9.0/2.0; - break; - case 2: a[0]=0.; a[1]=9.0; a[2]=-45.0/2.0; a[3]=27.0/2.0; break; - case 3: + case 2: a[0]=0.; a[1]=-9.0/2.0; a[2]=18.0; a[3]=-27.0/2.0; break; + case 3: + a[0]=0.; + a[1]=1.; + a[2]=-9.0/2.0; + a[3]=9.0/2.0; + break; default: Assert(false, ExcInternalError()); } @@ -142,33 +142,33 @@ vector LagrangeEquidistant::compute_coefficients(unsigned int n, unsigne a[4]=32.0/3.0; break; case 1: - a[0]=0.; - a[1]=-1.; - a[2]=22.0/3.0; - a[3]=-16.0; - a[4]=32.0/3.0; - break; - case 2: a[0]=0.; a[1]=16.0; a[2]=-208.0/3.0; a[3]=96.0; a[4]=-128.0/3.0; break; - case 3: + case 2: a[0]=0.; a[1]=-12.0; a[2]=76.0; a[3]=-128.0; a[4]=64.0; break; - case 4: + case 3: a[0]=0.; a[1]=16.0/3.0; a[2]=-112.0/3.0; a[3]=224.0/3.0; a[4]=-128.0/3.0; break; + case 4: + a[0]=0.; + a[1]=-1.; + a[2]=22.0/3.0; + a[3]=-16.0; + a[4]=32.0/3.0; + break; default: Assert(false, ExcInternalError()); } diff --git a/tests/base/polynomial_test.cc b/tests/base/polynomial_test.cc index 1a7af8f700..bdc4267e74 100644 --- a/tests/base/polynomial_test.cc +++ b/tests/base/polynomial_test.cc @@ -47,39 +47,23 @@ int main(int, char) LagrangeEquidistant polynom(order, s_point); // support points in vertices - if (order>0) - for (unsigned int i=0; i<=1; ++i) - { - double x=i; - polynom.value(x, values); - deallog << " p_" << s_point << "(" << x << ")"; -// deallog << "=" << values[0]; - if (equals_delta_ij(values[0], s_point, i)) - deallog << " ok"; - else - deallog << " false"; - deallog << endl; - } - // support points on line - if (order>1) - for (unsigned int i=1; i(i)/order; - polynom.value(x, values); - deallog << " p_" << s_point << "(" << x << ")"; -// deallog << "=" << values[0]; - if (equals_delta_ij(values[0], s_point, i+1)) - deallog << " ok"; - else - deallog << " false"; - deallog << endl; - } - deallog << endl; + for (unsigned int i=0; i<=order; ++i) + { + double x=static_cast(i)/order; + polynom.value(x, values); + deallog << " p_" << s_point << "(" << x << ")"; +// deallog << "=" << values[0]; + if (equals_delta_ij(values[0], s_point, i)) + deallog << " ok"; + else + deallog << " false"; + deallog << endl; + } } } deallog << endl << "Test derivatives computed by the Horner scheme:" << endl; - LagrangeEquidistant pol(4, 3); + LagrangeEquidistant pol(4, 2); vector v_horner(6); for (unsigned int i=0; i<=10; ++i) { @@ -87,8 +71,8 @@ int main(int, char) deallog << "x=" << xi << ", all derivatives: "; vector v_exact(6); - v_exact[0]=64.0*xi*xi*xi*xi-128.0*xi*xi*xi+76.0*xi*xi-12.0*xi; v_exact[1]=256.0*xi*xi*xi-384.0*xi*xi+152.0*xi-12.0; + v_exact[0]=64.0*xi*xi*xi*xi-128.0*xi*xi*xi+76.0*xi*xi-12.0*xi; v_exact[2]=768.0*xi*xi-768.0*xi+152.0; v_exact[3]=1536*xi-768; v_exact[4]=1536; diff --git a/tests/base/polynomial_test.checked b/tests/base/polynomial_test.checked index 06dbd848bb..71d3380875 100644 --- a/tests/base/polynomial_test.checked +++ b/tests/base/polynomial_test.checked @@ -3,75 +3,61 @@ DEAL::LagrangeEquidistant polynoms: DEAL::Polynomial p of order 1 DEAL:: p_0(0.00000) ok DEAL:: p_0(1.00000) ok - DEAL:: p_1(0.00000) ok DEAL:: p_1(1.00000) ok - DEAL::Polynomial p of order 2 DEAL:: p_0(0.00000) ok -DEAL:: p_0(1.00000) ok DEAL:: p_0(0.500000) ok - +DEAL:: p_0(1.00000) ok DEAL:: p_1(0.00000) ok -DEAL:: p_1(1.00000) ok DEAL:: p_1(0.500000) ok - +DEAL:: p_1(1.00000) ok DEAL:: p_2(0.00000) ok -DEAL:: p_2(1.00000) ok DEAL:: p_2(0.500000) ok - +DEAL:: p_2(1.00000) ok DEAL::Polynomial p of order 3 DEAL:: p_0(0.00000) ok -DEAL:: p_0(1.00000) ok DEAL:: p_0(0.333333) ok DEAL:: p_0(0.666667) ok - +DEAL:: p_0(1.00000) ok DEAL:: p_1(0.00000) ok -DEAL:: p_1(1.00000) ok DEAL:: p_1(0.333333) ok DEAL:: p_1(0.666667) ok - +DEAL:: p_1(1.00000) ok DEAL:: p_2(0.00000) ok -DEAL:: p_2(1.00000) ok DEAL:: p_2(0.333333) ok DEAL:: p_2(0.666667) ok - +DEAL:: p_2(1.00000) ok DEAL:: p_3(0.00000) ok -DEAL:: p_3(1.00000) ok DEAL:: p_3(0.333333) ok DEAL:: p_3(0.666667) ok - +DEAL:: p_3(1.00000) ok DEAL::Polynomial p of order 4 DEAL:: p_0(0.00000) ok -DEAL:: p_0(1.00000) ok DEAL:: p_0(0.250000) ok DEAL:: p_0(0.500000) ok DEAL:: p_0(0.750000) ok - +DEAL:: p_0(1.00000) ok DEAL:: p_1(0.00000) ok -DEAL:: p_1(1.00000) ok DEAL:: p_1(0.250000) ok DEAL:: p_1(0.500000) ok DEAL:: p_1(0.750000) ok - +DEAL:: p_1(1.00000) ok DEAL:: p_2(0.00000) ok -DEAL:: p_2(1.00000) ok DEAL:: p_2(0.250000) ok DEAL:: p_2(0.500000) ok DEAL:: p_2(0.750000) ok - +DEAL:: p_2(1.00000) ok DEAL:: p_3(0.00000) ok -DEAL:: p_3(1.00000) ok DEAL:: p_3(0.250000) ok DEAL:: p_3(0.500000) ok DEAL:: p_3(0.750000) ok - +DEAL:: p_3(1.00000) ok DEAL:: p_4(0.00000) ok -DEAL:: p_4(1.00000) ok DEAL:: p_4(0.250000) ok DEAL:: p_4(0.500000) ok DEAL:: p_4(0.750000) ok - +DEAL:: p_4(1.00000) ok DEAL::Test derivatives computed by the Horner scheme: DEAL::x=0.00000, all derivatives: ok -- 2.39.5