Line data Source code
1 0 : // Distributed under the MIT License.
2 : // See LICENSE.txt for details.
3 :
4 : #pragma once
5 :
6 : #include <cstddef>
7 :
8 : /// \cond
9 : class DataVector;
10 : class Matrix;
11 : namespace Spectral {
12 : enum class Quadrature : uint8_t;
13 : } // namespace Spectral
14 : /// \endcond
15 :
16 : namespace Spectral {
17 :
18 : /*!
19 : * \ingroup SpectralGroup
20 : *
21 : * \brief A collection of helper functions for Chebyshev polynomials
22 : *
23 : * \details The Chebyshev polynomials are given by:
24 : * \f[
25 : * T_n (x) = \cos(n \theta), \text{ } \theta = \arccos(x)
26 : * \f]
27 : *
28 : * The Chebyshev expansion of a function \f$f \in [-1,1] \f$ is given by:
29 : * \f[
30 : * f(x) = \sum_{n=0}^{\infty} f_n T_n(x)
31 : * \f]
32 : * where
33 : * \f[
34 : * f_n = \frac{1}{c_n} \int_{-1}^1 f(x) T_n(x) w(x) dx
35 : * \f]
36 : * where the weight function is \f$w(x) = (1-x^2)^{-1/2}\f$ and the basis
37 : * function normalization value is given by:
38 : * \f{align*}{
39 : * c_n
40 : * &=\begin{cases}
41 : * \hfil \pi \hfil & \text{if } n = 0 \\
42 : * \hfil \frac{\pi}{2} \hfil & \text{otherwise}
43 : * \end{cases}
44 : * \f}
45 : *
46 : * If a function is discretized at \f$N+1\f$ collocation points, the modal
47 : * representation will have \f$N+1\f$ spectral coefficients consisting of
48 : * \f[
49 : * f_n \qquad \text{for } n = 0, \ldots, N
50 : * \f]
51 : *
52 : * For more details about using Chebyshev polynomials see e.g. \cite Boyd2001
53 : * and \cite Fornberg1996.
54 : */
55 1 : class Chebyshev {
56 : public:
57 : /*!
58 : * \brief Value of the basis function \f$\Phi_k(x) = T_k(x)\f$
59 : */
60 : template <typename T>
61 1 : static T basis_function_value(size_t k, const T& x);
62 :
63 : /*!
64 : * \brief The normalization square \f$c_k\f$ of the basis function
65 : * \f$\Phi_k(x)\f$, i.e. the definite integral of its square.
66 : */
67 1 : static double basis_function_normalization_square(size_t k);
68 :
69 : /*!
70 : * \brief Collocation points \f$\{x_i\}\f$
71 : *
72 : * \details The collocation points on the interval \f$[-1, 1]\f$ are given by
73 : * \f{align*}{
74 : * x_i
75 : * &=\begin{cases}
76 : * \hfil - \cos \frac{(2i+1)\pi}{2N+2} \hfil &
77 : * \text{for Quadrature::Gauss} \\
78 : * \hfil - \cos \frac{i \pi}{N} \hfil &
79 : * \text{for Quadrature::GaussLobatto}
80 : * \end{cases}
81 : * \f}
82 : */
83 : template <Quadrature quadrature>
84 1 : static DataVector collocation_points(size_t num_points);
85 :
86 : /*!
87 : * \brief Integration weights \f$\{w_i\}\f$
88 : *
89 : * \details The integration weights are used to approximate the weighted
90 : * integral \f$Q[f]=\int f(x)w(x)\mathrm{d}x\approx \sum_k f_k w_k\f$.
91 : * For Quadrature::Gauss, the weights are given by:
92 : * \f[
93 : * w_i = \frac{\pi}{N+1}
94 : * \f]
95 : *
96 : * For Quadrature::GaussLobatto, the weights are given by:
97 : * \f{align*}{
98 : * w_i
99 : * &=\begin{cases}
100 : * \hfil \frac{\pi}{2N} \hfil & \text{for } j = 0, N\\
101 : * \hfil \frac{\pi}{N} \hfil & \text{for } j = 1, \ldots, N-1
102 : * \end{cases}
103 : * \f}
104 : *
105 : * \note These weights are used to compute the modes from the nodal
106 : * values. They are not used to evaluate definite or indefinite
107 : * integrals directly from the nodal values.
108 : */
109 : template <Quadrature quadrature>
110 1 : static DataVector integration_weights(size_t num_points);
111 :
112 : /*!
113 : * \brief Matrix used to compute the modes of the indefinite
114 : * integral from the modes of the integrand such that the constant
115 : * of integration is determined by requiring the integral to be zero
116 : * at \f$x=-1\f$.
117 : *
118 : * \details Chebyshev polynomials satisfy the identity:
119 : * \f[
120 : * \int^x dy T_n (y) = \left\{ \frac{T_{n+1}(x)}{2(n+1)} +
121 : * \frac{T_{n-1}(x)}{2(n-1)} \right\}, \qquad n \geq 2
122 : * \f]
123 : *
124 : * Thus the modes \f$\tilde{f}_j\f$ of the integral are given as:
125 : * \f{align*}{
126 : * \tilde{f}_i
127 : * &=\begin{cases}
128 : * \hfil \frac{f_{i-1}-f_{i+1}}{2i}, \hfil & \text{for } i > 1\\
129 : * \hfil f_0 - \frac{f_2}{2} \hfil & \text{for } i = 1
130 : * \end{cases}
131 : * \f}
132 : * where \f$f_j\f$ are the modes of the integrand.
133 : *\f$\tilde{f}_0\f$ is a constant of integration, which we choose such that
134 : * the integral is 0 at the left boundary of the domain (\f$x=-1\f$). The
135 : * condition for this is:
136 : *
137 : * \f[
138 : * \tilde{f}_0=\sum_{i=1}^{N}(-1)^{i+1}\tilde{f}_i
139 : * \f]
140 : */
141 1 : static Matrix indefinite_integral_matrix(size_t num_points);
142 :
143 : /*!
144 : * \brief Row-vector used to compute the definite integral from the modes of
145 : * the integrand
146 : *
147 : * \details Given the modes \f$\tilde{f}_j\f$ of the integral (see
148 : * indefinite_integral_matrix), the definite integral is given by evaluating
149 : * the series at \f$x=1\f$ and \f$x=-1\f$ and taking the difference. Given
150 : * that \f$T_n(1) = 1\f$ and \f$T_n(-1) = (-1)^n\f$, this means multiplying
151 : * the indefinite_integral_matrix by the row-vector \f$\{0, 2, 0, 2, \ldots
152 : * \}\f$ which yields:
153 : * \f{align*}{
154 : * \tilde{q}_i
155 : * &=\begin{cases}
156 : * \hfil -\frac{2}{n^2-1}, \hfil & \text{for } i \text{ even}\\
157 : * \hfil 0 \hfil & \text{for } i \text{ odd}
158 : * \end{cases}
159 : * \f}
160 : */
161 1 : static Matrix definite_integral_matrix(size_t num_points);
162 : };
163 : } // namespace Spectral
|
