#include <vector>
+template <int dim> class Point;
/**
* A namespace in which classes relating to the description of
*/
Polynomial<number>& operator *= (const double s);
+ /**
+ * Multiply with another polynomial.
+ */
+ Polynomial<number>& operator *= (const Polynomial<number>& p);
+
/**
* Add a second polynomial.
*/
const unsigned int support_point);
};
-
+/**
+ * Lagrange polynomials for an arbistrary set of interpolation points.
+ *
+ * @author Guido Kanschat, 2005
+ */
+ class Lagrange
+ {
+ public:
+ /**
+ * Given a set of points, this
+ * function returns all
+ * Lagrange polynomials for
+ * interpolation of these
+ * points. The number of
+ * polynomials is equal to the
+ * number of points and the
+ * maximum degree is one less.
+ */
+ static
+ std::vector<Polynomial<double> >
+ generate_complete_basis (const std::vector<Point<1> >& points);
+ };
+
+
+
/**
* Legendre polynomials of arbitrary degree on <tt>[0,1]</tt>.
*
#include <base/polynomial.h>
+#include <base/point.h>
#include <base/exceptions.h>
#include <base/thread_management.h>
}
+ template <typename number>
+ Polynomial<number>&
+ Polynomial<number>::operator *= (const Polynomial<number>& p)
+ {
+ // Degree of the product
+ unsigned int new_degree = this->degree() + p.degree();
+
+ std::vector<number> new_coefficients(new_degree+1, 0.);
+
+ for (unsigned int i=0; i<p.coefficients.size(); ++i)
+ for (unsigned int j=0; j<this->coefficients.size(); ++j)
+ new_coefficients[i+j] += this->coefficients[j]*p.coefficients[i];
+ this->coefficients = new_coefficients;
+
+ return *this;
+ }
+
+
template <typename number>
Polynomial<number>&
Polynomial<number>::operator += (const Polynomial<number>& p)
}
+//----------------------------------------------------------------------//
+
+
+ std::vector<Polynomial<double> >
+ Lagrange::generate_complete_basis (const std::vector<Point<1> >& points)
+ {
+ std::vector<Polynomial<double> > p(points.size());
+ // polynomials are built as
+ // products of linear
+ // factors. The coefficient in
+ // front of the linear term is
+ // always 1.
+ std::vector<double> linear(2, 1.);
+ // We start with a constant polynomial
+ std::vector<double> one(1, 1.);
+
+ for (unsigned int i=0;i<p.size();++i)
+ {
+ // Construct interpolation formula
+ p[i] = Polynomial<double>(one);
+ for (unsigned int k=0;k<p.size();++k)
+ if (k != i)
+ {
+ linear[0] = -points[k](0);
+ Polynomial<double> factor(linear);
+ factor *= 1./(points[i](0)-points[k](0));
+ p[i] *= factor;
+ }
+ }
+ return p;
+ }
+
// ------------------ class Legendre --------------- //