]> https://gitweb.dealii.org/ - dealii-svn.git/commitdiff
Lagrange interpolation for arbitrary points
authorguido <guido@0785d39b-7218-0410-832d-ea1e28bc413d>
Mon, 30 May 2005 14:46:39 +0000 (14:46 +0000)
committerguido <guido@0785d39b-7218-0410-832d-ea1e28bc413d>
Mon, 30 May 2005 14:46:39 +0000 (14:46 +0000)
git-svn-id: https://svn.dealii.org/trunk@10778 0785d39b-7218-0410-832d-ea1e28bc413d

deal.II/base/source/polynomial.cc
tests/base/polynomial1d.cc
tests/base/polynomial_lagrange.cc [new file with mode: 0644]
tests/results/i686-pc-linux-gnu+gcc3.2/base/polynomial_lagrange.output [new file with mode: 0644]

index 1e3e6589bded07df0afb14781066ed3e4c7612e2..b57c923841aedc1aa3b7df9bf91dbaa557fa3ade 100644 (file)
@@ -628,7 +628,7 @@ namespace Polynomials
       {
                                         // Construct interpolation formula
        p[i] = Polynomial<double>(one);
-       for (unsigned int k=0;k<p.size();++k)
+       for (unsigned int k=0;k<points.size();++k)
          if (k != i)
            {
              linear[0] = -points[k](0);
@@ -637,6 +637,30 @@ namespace Polynomials
              p[i] *= factor;
            }
       }
+
+                                    // Since the previous operation
+                                    // is subject to round-off error
+                                    // amplification, we correct the
+                                    // polynomials here.
+    for (unsigned int i=0;i<p.size();++i)
+      {
+       for (unsigned int k=0;k<points.size();++k)
+         {
+           const double value = p[i].value(points[k](0));
+           Polynomial<double> q = p[k];
+           if (i==k)
+             {
+               q *= 1.-value;
+               p[i] += q;
+             }
+           else
+             {
+               q *= -value;
+               p[i] += q;
+             }
+         }
+      }
+    
     return p;
   }
   
index e2a5f782508c551065ad3e46f7404928088c1d4e..9f59c9778356beb1cb74ecf80496a3d8043e3c6a 100644 (file)
@@ -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<p.size();++i)
     for (unsigned int j=0;j<p.size();++j)
diff --git a/tests/base/polynomial_lagrange.cc b/tests/base/polynomial_lagrange.cc
new file mode 100644 (file)
index 0000000..ccea302
--- /dev/null
@@ -0,0 +1,128 @@
+//-----------------------------------------------------------------------------
+//    $Id$
+//    Version: $Name$ 
+//
+//    Copyright (C) 2000, 2001, 2002, 2003, 2004, 2005 by the deal.II authors
+//
+//    This file is subject to QPL and may not be  distributed
+//    without copyright and license information. Please refer
+//    to the file deal.II/doc/license.html for the  text  and
+//    further information on this license.
+//
+//-----------------------------------------------------------------------------
+
+// Test Lagrange interpolation
+
+#include "../tests.h"
+#include <iostream>
+#include <fstream>
+#include <cmath>
+
+#include <base/logstream.h>
+#include <base/polynomial.h>
+#include <base/quadrature_lib.h>
+
+using namespace Polynomials;
+
+void check_interpolation (const std::vector<Polynomial<double> >& p,
+                         const std::vector<Point<1> >& x)
+{
+  for (unsigned int i=0;i<p.size();++i)
+    {
+      deallog << i;
+      for (unsigned int k=0;k<x.size();++k)
+       {
+         deallog << '.';
+         const double y = p[i].value(x[k](0));
+         if (i == k)
+           {
+             if (std::fabs(y-1.) > 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<Polynomial<double> > 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<Polynomial<double> > p = LagrangeEquidistant::generate_complete_basis(n);
+  std::vector<Point<1> > 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 (file)
index 0000000..0db05e9
--- /dev/null
@@ -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..........

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.