From 4b50d27eb3372569dc42b6eecdc5bfdeac5d3a00 Mon Sep 17 00:00:00 2001
From: guido <guido@0785d39b-7218-0410-832d-ea1e28bc413d>
Date: Mon, 30 May 2005 14:46:39 +0000
Subject: [PATCH] Lagrange interpolation for arbitrary points

git-svn-id: https://svn.dealii.org/trunk@10778 0785d39b-7218-0410-832d-ea1e28bc413d
---
 deal.II/base/source/polynomial.cc             |  26 +++-
 tests/base/polynomial1d.cc                    |   2 +-
 tests/base/polynomial_lagrange.cc             | 128 ++++++++++++++++++
 .../base/polynomial_lagrange.output           | 116 ++++++++++++++++
 4 files changed, 270 insertions(+), 2 deletions(-)
 create mode 100644 tests/base/polynomial_lagrange.cc
 create mode 100644 tests/results/i686-pc-linux-gnu+gcc3.2/base/polynomial_lagrange.output

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<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;
   }
   
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<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
index 0000000000..ccea302fbf
--- /dev/null
+++ b/tests/base/polynomial_lagrange.cc
@@ -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
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..........
-- 
2.39.5