Reference documentation for deal.II version 9.0.0
polynomials_integrated_legendre_sz.cc
1 // ---------------------------------------------------------------------
2 //
3 // Copyright (C) 2015 - 2017 by the deal.II authors
4 //
5 // This file is part of the deal.II library.
6 //
7 // The deal.II library is free software; you can use it, redistribute
8 // it, and/or modify it under the terms of the GNU Lesser General
9 // Public License as published by the Free Software Foundation; either
10 // version 2.1 of the License, or (at your option) any later version.
11 // The full text of the license can be found in the file LICENSE at
12 // the top level of the deal.II distribution.
13 //
14 // ---------------------------------------------------------------------
15 
16 
17 #include <deal.II/base/polynomials_integrated_legendre_sz.h>
18 
19 DEAL_II_NAMESPACE_OPEN
20 
22  :
23  Polynomials::Polynomial<double> (get_coefficients(k))
24 {}
25 
26 
27 
28 const std::vector<double> IntegratedLegendreSZ::get_coefficients (const unsigned int k)
29 {
30  std::vector<double> coefficients(k+1);
31 
32  // first two polynomials are hard-coded:
33  if (k==0)
34  {
35  coefficients[0] = -1.;
36  return coefficients;
37  }
38  else if (k==1)
39  {
40  coefficients[0] = 0.;
41  coefficients[1] = 1.;
42  return coefficients;
43  }
44 
45  // General formula is:
46  // k*L_{k}(x) = (2*k-3)*x*L_{k-1} - (k-3)*L_{k-2}.
47  std::vector<double> coefficients_km2 = get_coefficients(k-2);
48  std::vector<double> coefficients_km1 = get_coefficients(k-1);
49 
50  const double a = 1.0 / k;
51  const double b = 2.0*k - 3.0;
52  const double c = k - 3.0;
53 
54  // To maintain stability, delay the division (multiplication by a) until the end.
55  for (unsigned int i=1; i<=k-2; i++)
56  {
57  coefficients[i] = b*coefficients_km1[i-1] - c*coefficients_km2[i];
58  }
59 
60  coefficients[0] = -c*coefficients_km2[0];
61  coefficients[k] = b*coefficients_km1[k-1];
62  coefficients[k-1] = b*coefficients_km1[k-2];
63 
64  for (unsigned int i=0; i<coefficients.size(); i++)
65  {
66  coefficients[i] *= a;
67  }
68 
69  return coefficients;
70 }
71 
72 
73 
74 std::vector<Polynomials::Polynomial<double> >
76 {
77  std::vector<Polynomials::Polynomial<double> > v;
78  v.reserve(degree + 1);
79  for (unsigned int i=0; i<=degree; ++i)
80  {
81  v.push_back (IntegratedLegendreSZ(i));
82  }
83  return v;
84 }
85 
86 
87 
88 DEAL_II_NAMESPACE_CLOSE
static std::vector< Polynomials::Polynomial< double > > generate_complete_basis(const unsigned int degree)
unsigned int degree() const
Definition: polynomial.h:769
static const std::vector< double > get_coefficients(const unsigned int k)
std::vector< double > coefficients
Definition: polynomial.h:250