Loading [MathJax]/extensions/TeX/newcommand.js
 deal.II version GIT relicensing-3124-g7c413cc06e 2025-04-23 12:20:01+00:00
\newcommand{\dealvcentcolon}{\mathrel{\mathop{:}}} \newcommand{\dealcoloneq}{\dealvcentcolon\mathrel{\mkern-1.2mu}=} \newcommand{\jump}[1]{\left[\!\left[ #1 \right]\!\right]} \newcommand{\average}[1]{\left\{\!\left\{ #1 \right\}\!\right\}}
All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Modules Pages Concepts
polynomials_integrated_legendre_sz.cc
Go to the documentation of this file.
1// ------------------------------------------------------------------------
2//
3// SPDX-License-Identifier: LGPL-2.1-or-later
4// Copyright (C) 2017 - 2023 by the deal.II authors
5//
6// This file is part of the deal.II library.
7//
8// Part of the source code is dual licensed under Apache-2.0 WITH
9// LLVM-exception OR LGPL-2.1-or-later. Detailed license information
10// governing the source code and code contributions can be found in
11// LICENSE.md and CONTRIBUTING.md at the top level directory of deal.II.
12//
13// ------------------------------------------------------------------------
14
15
17
18#include <vector>
19
21
23 : Polynomials::Polynomial<double>(get_coefficients(k))
24{}
25
26
27
28std::vector<double>
30{
31 std::vector<double> coefficients(k + 1);
32
33 // first two polynomials are hard-coded:
34 if (k == 0)
35 {
36 coefficients[0] = -1.;
37 return coefficients;
38 }
39 else if (k == 1)
40 {
41 coefficients[0] = 0.;
42 coefficients[1] = 1.;
43 return coefficients;
44 }
45
46 // General formula is:
47 // k*L_{k}(x) = (2*k-3)*x*L_{k-1} - (k-3)*L_{k-2}.
48 std::vector<double> coefficients_km2 = get_coefficients(k - 2);
49 std::vector<double> coefficients_km1 = get_coefficients(k - 1);
50
51 const double a = 1.0 / k;
52 const double b = 2.0 * k - 3.0;
53 const double c = k - 3.0;
54
55 // To maintain stability, delay the division (multiplication by a) until the
56 // end.
57 for (unsigned int i = 1; i <= k - 2; ++i)
58 {
59 coefficients[i] = b * coefficients_km1[i - 1] - c * coefficients_km2[i];
60 }
61
62 coefficients[0] = -c * coefficients_km2[0];
63 coefficients[k] = b * coefficients_km1[k - 1];
64 coefficients[k - 1] = b * coefficients_km1[k - 2];
65
66 for (double &coefficient : coefficients)
67 {
68 coefficient *= a;
69 }
70
71 return coefficients;
72}
73
74
75
76std::vector<Polynomials::Polynomial<double>>
78{
79 std::vector<Polynomials::Polynomial<double>> v;
80 v.reserve(degree + 1);
81 for (unsigned int i = 0; i <= degree; ++i)
82 {
83 v.push_back(IntegratedLegendreSZ(i));
84 }
85 return v;
86}
87
88
89
static std::vector< Polynomials::Polynomial< double > > generate_complete_basis(const unsigned int degree)
static std::vector< double > get_coefficients(const unsigned int k)
std::vector< double > coefficients
Definition polynomial.h:307
unsigned int degree() const
Definition polynomial.h:811
#define DEAL_II_NAMESPACE_OPEN
Definition config.h:35
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:36