From: Guido Kanschat Date: Mon, 30 May 2005 14:46:39 +0000 (+0000) Subject: Lagrange interpolation for arbitrary points X-Git-Tag: v8.0.0~13778 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=85a094b9548e8a7b9a31dca57ff0b0ff14b62c0d;p=dealii.git Lagrange interpolation for arbitrary points git-svn-id: https://svn.dealii.org/trunk@10778 0785d39b-7218-0410-832d-ea1e28bc413d --- diff --git a/deal.II/base/source/polynomial.cc b/deal.II/base/source/polynomial.cc index 1e3e6589bd..b57c923841 100644 --- a/deal.II/base/source/polynomial.cc +++ b/deal.II/base/source/polynomial.cc @@ -628,7 +628,7 @@ namespace Polynomials { // Construct interpolation formula p[i] = Polynomial(one); - for (unsigned int k=0;k q = p[k]; + if (i==k) + { + q *= 1.-value; + p[i] += q; + } + else + { + q *= -value; + p[i] += q; + } + } + } + return p; } diff --git a/tests/base/polynomial1d.cc b/tests/base/polynomial1d.cc index e2a5f78250..9f59c97783 100644 --- a/tests/base/polynomial1d.cc +++ b/tests/base/polynomial1d.cc @@ -127,7 +127,7 @@ int main () q.push_back(LagrangeEquidistant(6, i)); } - // We add 1.0001 bacuse of bugs in + // We add 1.0001 because of bugs in // the ostream classes for (unsigned int i=0;i +#include +#include + +#include +#include +#include + +using namespace Polynomials; + +void check_interpolation (const std::vector >& p, + const std::vector >& x) +{ + for (unsigned int i=0;i 2.e-10) + deallog << "Error1 lg y=" << std::log10(std::fabs(y-1.)) + << std::endl; + } + else + { + if (std::fabs(y) > 2.e-10) + deallog << "Error0 lg y=" << std::log10(std::fabs(y)) + << std::endl; + } + } + deallog << std::endl; + } +} + + +void +check_poly (const Quadrature<1>& q) +{ + deallog << "Points: " << q.n_quadrature_points << std::endl; + std::vector > p = Lagrange::generate_complete_basis(q.get_points()); + check_interpolation(p, q.get_points()); +} + + +void +check_lge (unsigned int n) +{ + deallog << "Points: " << n+1 << std::endl; + std::vector > p = LagrangeEquidistant::generate_complete_basis(n); + std::vector > x(n+1); + const double h = 1./n; + for (unsigned int i=0;i<=n;++i) + x[i](0) = h*i; + check_interpolation(p, x); +} + + +int main() +{ + std::ofstream logfile("polynomial_lagrange.output"); + logfile.precision(3); + deallog.attach(logfile); + deallog.depth_console(0); +// deallog.threshold_double(1.e-10); + + QTrapez<1> trapez; + QSimpson<1> simpson; + QIterated<1> equi7(trapez, 6); + QIterated<1> equi10(trapez, 9); + + QGauss<1> g2(2); + QGauss<1> g3(3); + QGauss<1> g7(7); + QGauss<1> g10(10); + + QGaussLobatto<1> gl2(2); + QGaussLobatto<1> gl3(3); + QGaussLobatto<1> gl4(4); + QGaussLobatto<1> gl5(5); + QGaussLobatto<1> gl6(6); + QGaussLobatto<1> gl7(7); + QGaussLobatto<1> gl10(10); + + deallog.push("LagrangeEquidistant"); + check_lge(6); + check_lge(9); + deallog.pop(); + deallog.push("Equidistant"); + check_poly(trapez); + check_poly(simpson); + check_poly(equi7); + check_poly(equi10); + deallog.pop(); + deallog.push("Gauss"); + check_poly(g2); + check_poly(g3); + check_poly(g7); + check_poly(g10); + deallog.pop(); + deallog.push("GaussLobatto"); + check_poly(gl2); + check_poly(gl3); + check_poly(gl4); + check_poly(gl5); + check_poly(gl6); + check_poly(gl7); + check_poly(gl10); + deallog.pop(); +}