From: guido Date: Mon, 30 May 2005 14:46:39 +0000 (+0000) Subject: Lagrange interpolation for arbitrary points X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=4b50d27eb3372569dc42b6eecdc5bfdeac5d3a00;p=dealii-svn.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(); +} diff --git a/tests/results/i686-pc-linux-gnu+gcc3.2/base/polynomial_lagrange.output b/tests/results/i686-pc-linux-gnu+gcc3.2/base/polynomial_lagrange.output new file mode 100644 index 0000000000..0db05e9f74 --- /dev/null +++ b/tests/results/i686-pc-linux-gnu+gcc3.2/base/polynomial_lagrange.output @@ -0,0 +1,116 @@ + +DEAL:LagrangeEquidistant::Points: 7 +DEAL:LagrangeEquidistant::0....... +DEAL:LagrangeEquidistant::1....... +DEAL:LagrangeEquidistant::2....... +DEAL:LagrangeEquidistant::3....... +DEAL:LagrangeEquidistant::4....... +DEAL:LagrangeEquidistant::5....... +DEAL:LagrangeEquidistant::6....... +DEAL:LagrangeEquidistant::Points: 10 +DEAL:LagrangeEquidistant::0.......... +DEAL:LagrangeEquidistant::1.......... +DEAL:LagrangeEquidistant::2.......... +DEAL:LagrangeEquidistant::3.......... +DEAL:LagrangeEquidistant::4.......... +DEAL:LagrangeEquidistant::5.......... +DEAL:LagrangeEquidistant::6.......... +DEAL:LagrangeEquidistant::7.......... +DEAL:LagrangeEquidistant::8.......... +DEAL:LagrangeEquidistant::9.......... +DEAL:Equidistant::Points: 2 +DEAL:Equidistant::0.. +DEAL:Equidistant::1.. +DEAL:Equidistant::Points: 3 +DEAL:Equidistant::0... +DEAL:Equidistant::1... +DEAL:Equidistant::2... +DEAL:Equidistant::Points: 7 +DEAL:Equidistant::0....... +DEAL:Equidistant::1....... +DEAL:Equidistant::2....... +DEAL:Equidistant::3....... +DEAL:Equidistant::4....... +DEAL:Equidistant::5....... +DEAL:Equidistant::6....... +DEAL:Equidistant::Points: 10 +DEAL:Equidistant::0.......... +DEAL:Equidistant::1.......... +DEAL:Equidistant::2.......... +DEAL:Equidistant::3.......... +DEAL:Equidistant::4.......... +DEAL:Equidistant::5.......... +DEAL:Equidistant::6.......... +DEAL:Equidistant::7.......... +DEAL:Equidistant::8.......... +DEAL:Equidistant::9.......... +DEAL:Gauss::Points: 2 +DEAL:Gauss::0.. +DEAL:Gauss::1.. +DEAL:Gauss::Points: 3 +DEAL:Gauss::0... +DEAL:Gauss::1... +DEAL:Gauss::2... +DEAL:Gauss::Points: 7 +DEAL:Gauss::0....... +DEAL:Gauss::1....... +DEAL:Gauss::2....... +DEAL:Gauss::3....... +DEAL:Gauss::4....... +DEAL:Gauss::5....... +DEAL:Gauss::6....... +DEAL:Gauss::Points: 10 +DEAL:Gauss::0.......... +DEAL:Gauss::1.......... +DEAL:Gauss::2.......... +DEAL:Gauss::3.......... +DEAL:Gauss::4.......... +DEAL:Gauss::5.......... +DEAL:Gauss::6.......... +DEAL:Gauss::7.......... +DEAL:Gauss::8.......... +DEAL:Gauss::9.......... +DEAL:GaussLobatto::Points: 2 +DEAL:GaussLobatto::0.. +DEAL:GaussLobatto::1.. +DEAL:GaussLobatto::Points: 3 +DEAL:GaussLobatto::0... +DEAL:GaussLobatto::1... +DEAL:GaussLobatto::2... +DEAL:GaussLobatto::Points: 4 +DEAL:GaussLobatto::0.... +DEAL:GaussLobatto::1.... +DEAL:GaussLobatto::2.... +DEAL:GaussLobatto::3.... +DEAL:GaussLobatto::Points: 5 +DEAL:GaussLobatto::0..... +DEAL:GaussLobatto::1..... +DEAL:GaussLobatto::2..... +DEAL:GaussLobatto::3..... +DEAL:GaussLobatto::4..... +DEAL:GaussLobatto::Points: 6 +DEAL:GaussLobatto::0...... +DEAL:GaussLobatto::1...... +DEAL:GaussLobatto::2...... +DEAL:GaussLobatto::3...... +DEAL:GaussLobatto::4...... +DEAL:GaussLobatto::5...... +DEAL:GaussLobatto::Points: 7 +DEAL:GaussLobatto::0....... +DEAL:GaussLobatto::1....... +DEAL:GaussLobatto::2....... +DEAL:GaussLobatto::3....... +DEAL:GaussLobatto::4....... +DEAL:GaussLobatto::5....... +DEAL:GaussLobatto::6....... +DEAL:GaussLobatto::Points: 10 +DEAL:GaussLobatto::0.......... +DEAL:GaussLobatto::1.......... +DEAL:GaussLobatto::2.......... +DEAL:GaussLobatto::3.......... +DEAL:GaussLobatto::4.......... +DEAL:GaussLobatto::5.......... +DEAL:GaussLobatto::6.......... +DEAL:GaussLobatto::7.......... +DEAL:GaussLobatto::8.......... +DEAL:GaussLobatto::9..........