Reference documentation for deal.II version 9.1.1
\(\newcommand{\dealcoloneq}{\mathrel{\vcenter{:}}=}\)
polynomials_integrated_legendre_sz.cc
1 // ---------------------------------------------------------------------
2 //
3 // Copyright (C) 2015 - 2019 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.md at
12 // the top level directory of deal.II.
13 //
14 // ---------------------------------------------------------------------
15 
16 
17 #include <deal.II/base/polynomials_integrated_legendre_sz.h>
18 
19 DEAL_II_NAMESPACE_OPEN
20 
22  : Polynomials::Polynomial<double>(get_coefficients(k))
23 {}
24 
25 
26 
27 const std::vector<double>
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
55  // end.
56  for (unsigned int i = 1; i <= k - 2; i++)
57  {
58  coefficients[i] = b * coefficients_km1[i - 1] - c * coefficients_km2[i];
59  }
60 
61  coefficients[0] = -c * coefficients_km2[0];
62  coefficients[k] = b * coefficients_km1[k - 1];
63  coefficients[k - 1] = b * coefficients_km1[k - 2];
64 
65  for (double &coefficient : coefficients)
66  {
67  coefficient *= a;
68  }
69 
70  return coefficients;
71 }
72 
73 
74 
75 std::vector<Polynomials::Polynomial<double>>
77 {
78  std::vector<Polynomials::Polynomial<double>> v;
79  v.reserve(degree + 1);
80  for (unsigned int i = 0; i <= degree; ++i)
81  {
82  v.push_back(IntegratedLegendreSZ(i));
83  }
84  return v;
85 }
86 
87 
88 
89 DEAL_II_NAMESPACE_CLOSE
static std::vector< Polynomials::Polynomial< double > > generate_complete_basis(const unsigned int degree)
unsigned int degree() const
Definition: polynomial.h:775
static const std::vector< double > get_coefficients(const unsigned int k)
std::vector< double > coefficients
Definition: polynomial.h:264